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NOMENCLATURE 


English 

Symbol 

Definition 

Units Used 

a 

Propellant burning rate coefficient. 

in/sec-psi 11 

c 

Specific heat. 

in-lbf/lbm*F 

C 

V 

Coefficient of variation; i.e., the ratio of the 
standard deviation to the mean. 

— 

e ij 

Strain tensor. 

— 

E 

Modulus of elasticity. 

lbf/in 2 

F 

Thrust. 

lbf 

K 

Statistical confidence coefficient. 

— 

X t 

Total impulse. 

lbf-sec 

£ 

Length of propellant grain. 

in. 

m 

Mass or ratio of inside to outside radius. 

slugs or 

”G 

Mass of propellant gases generated per unit time. 

slugs/sec 

n 

Burning rate exponent or number of observations 
of a statistically distributed variable. 

— 

P 

Pressure . 

lbf/in 2 

r 

Burning rate. 

in/sec 

r 

c 

Correlation coefficient. 

— 

R 

Radius . 

in. 

s 

Standard deviation of a sample of a statistically 
distributed variable. 

units vary 

s 

o 

Square root of the second moment of a sample 
distribution about an assumed zero mean. 

units vary 

S 

Burning perimeter. 

in. 

t 

Time. 

sec. 
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Symbol 

Definition 

Units Used 

T 

Temperature in thermoelastic analysis. 

°R 

T Gr’ T ref 

Grain temperature and standard grain temperature, 
respectively. 

°F 

X 

Radial coordinate in thermoelastic analysis. 

in. 

X 

Value of general statistically distributed 
variable. 

units vary 

y.y’ 

Distance propellant has burned from initial 
lateral surface of circular perforated grain 
and from other initial surfaces, respectively. 

in. 

Greek 

Symbol 



a 

Linear coefficient of thermal expansion 

/°F 

3 

Volumetric coefficient of thermal expansion 

/°F 

*« 

The Kronecker delta (1 when i=j, 0 when i^j) 

— 

A 

Change or difference in quantity. 

units vary 

e 

Strain. 


n 

Compressibility . 

in 2 /lbf 

X 

Thermal conductivity. 

in-lbf/in-sec 

V 

Poisson's ratio. 

— 

P 

Density. 

slugs/in 3 

a 

The standard deviation of a statistically 
distributed variable; i.e., the square root of 
the second moment about its mean value. 

units vary 

0 

P 

Temperature sensitivity of burning rate at 
constant pressure . 

/°R 

T 

Thickness . 

in. 
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Subscripts 

av 

c 

e 

i 

max 

o 

P 

r 

s 

ss 

t 

w 

y 

z 

0 


Definition 

Average or mean value. 

Motor case or compressed state of propellant just outside 
the heat-affected zone. 

Grain exterior surface position. 

Grain exterior surface position or initial condition. 
Maximum value. 

Unpressurized and unheated state of propellant or in the 
thermoelastic analysis the ambient transient condition. 

Propellant. 

Radial coordinate. 

Surface condition. 

Steady state 

Transient 

Propellant web. 

Distance burned. 

Axial coordinate. 

Tangential coordinate. 


x 



I. INTRODUCTION AND SUMMARY 


This report presents the results of research performed at Auburn 
University during the period October 28, 1975, to September 30, 1976, 
under Modifications Nos. 17 and 18 to the Cooperative Agreement, dated 
February 11, 1969, between NASA Marshall Space Flight Center (MSFC) 
and Auburn University. The principal objective of the research was to 
further develop techniques for theoretical assessment of solid rocket 
motor (3RM) internal ballistic performance to include statistical in- 
vestigation of thrust imbalance of pairs of SRMs firing in parallel as 
on the booster stage of the Space Shuttle. 

The theoretical thrust imbalance of motor pairs has been previously 
investigated statistically by application of the Monte Carlo technique 
(Refs. 1-3). The results of this investigaticr. include a computer pro- 
gram which selects sets of significant variables on a probability basis 
and calculates the performance cnaracteristics for a large number of 
motor pairs using a mathematical model of the internal ballistics. Com- 
parison of such a statistical analysis for TITAN IIIC motor pairs with 
test results shows the theory underpredicts the standard deviation in 
maximum thrust imbalance by 20% with variability in burning times matched 
within 2% (Ref. 1). 

The referenced research disclosed a number of ways in which both 
nominal and off-nominal performance predictions for single and pairs of 
SRMs could be improved. The present report deals with the results of 
research to produce these improvements. 

Recent efforts by MSFC to describe thrust imbalance characteristics 
to be anticipated in a concise way that would be meaningful for systems 
performance studies have indicated that the Monte Carl ' program would be 
useful for this purpose. Two supplementary computer programs are included 
in this report to permit the critical imbalance parameters to be determined 
on a statistical basis throughout the operating times of the SRMs. The 
principal critical imbalance parameters are the thrust imbalance and its 
first time derivative. Data for these and other parameters of interest; 
e.g., impulse imbalance; arc generated by the Monte Carlo program and 
stored on magnetic, tape. One program is used to determine the maximum 
values of the parameters of interest during specified time periods. The 
second program computes statistical tolerance limits as a function of 
time based on the Monte Carlo sample. CalComp plots of the parameters 
of interest versus time are obtainable for both individual SRM pairs over- 
laid on each other and for the tolerance limits. These are illustrated 
in the report. Comparison of program results with flight test results 
for Titan IIIC are given and show good agreement. 

The correlation coefficients between thrust imbalance and its time 
derivative were determined at various times. The results show consider- 
able variation in the goodness of the correlation at various operating 
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times. A linear regression analysis shows little generalisation is 
possible as to the relation between these variables at the various 
operating times. On the other hand, there appears to be excellent cor- 
relation between thrust imbalance and the nominal slope of the thrust 
time trace which suggests several alternative approaches to that used in 
Refs. 1 and 2 for determining thrust imbalance. 

Another portion of this research treats certain effects of grain de- 
formation on the internal ballistics of SKMs . The report explains how the 
deformation can affect the burning rate of the propellant and account for 
a significant portion of the "scale factor" between large and small motor 
burning rate. The important factor is the tangential strain of the burn- 
ing surface of the propellant grain. A method developed by Smith (Ref. 4) 
has been used to estimate the strain at various times. The method was 
coupled with the design analysis computer program (Ref. 2) to analyze the 
effects on internal ballistics for circular-perforated grains. Comparisons 
are presented of the internal ballistic results with and without grain 
deformation and with actual test data for two different SRMs. The modified 
design analysis computer program is listed in the report along with a sample 
solution. Variations in the propellant and motor case properties (especially 
variation in the propellant modulus) could produce variations in the grain 
deformation causing each SRM of a pair to perform differently. However, 
further verification of the deformation analysis is in order before the 
results are incorporated into the Monte Carlo program. 

A third general aspect of the present investigation is that of changes 
in grain temperatures produced by high strain rates. Previous research 
(Ref. 2) indicated that thermoelastic coupling could have a significant 
effect upon the grain temperature distribution near the burning surface 
during highly transient chamber pressure conditions. This evaluation was 
based upon the assumption of a specified fixed surface temperature. How- 
ever, during transient conditions the surface temperature is not fixed 
but depends upon both the combustion chamber pressure and the thermo- 
elastic coupling itself which will alter the heat transfer from the com- 
bustion zone. Also, the previous analysis considered the effect of 
surface regression on the energy balance only in an approximate and in- 
tuitive way. In the present work, the problem was solved more rigorously. 

The surface regression term is included in the energy equation along with 
the strain rate (thermoolast ic) term. ’Ihc surface temperature is cal- 
culated by coupling the energy equation with a model of the combustion 
zone. This is accomplished using the Zeldovi ch-Novozhilov (Z-N) technique 
(Ref. 5) which relates the transient heat flux at the surface to steady 
state burning rate data. The Z-N method has the advantage over other 
approaches considered of not requiring detailed knowledge of the flame 
structure which would involve quantities that are difficult to measure. 

Solution of the differential equation for the energy balance of the 
solid propellant grain subject to the appropriate boundary conditions 
yields the temperature distribution and the surface regression rate. The 
solutions are obtained using, the numerical technique of Foster (Refs. 2 
and 6) applied to a circular-perforated case-bonded grain of finite 
length subjected to an increasing chamber pressure. 
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Comparisons are made of the results with and without thermoelastic 
coupling. It is concluded from the comparison that the effect of thermo- 
elastic coupling is very slight. However, it may not be insignificant. 
Additional research is indicated, particularly with respect to possible 
viscoelastic effects. 

Another facet of internal ballistic performance variation considered 
is the ability to assess the quality of ballistic predictions. It is im- 
portant to know the confidence with which the performance of an SRM as a 
function of time may be predicted at any phase cf a motor development 
program. In general, predictability will improve as more motors of a single 
type are manufactured and tested. As shown in this report, the Monte Carlo 
program provides an approach to assessing the degree of predictability 
both for single SRMs and pairs at various phases in the development pro- 
gram. However, further investigation has shown that it may prove impractical 
to obtain the type of design and manufacturing data necessary to establish 
predictability by this method. 

The final portion of the report identifies additions and other changes 
that have been made to the two computer programs listed in Ref. 2: the 

Monte Carlo program and the design analysis program. Each change is 
identified by page number and line, as listed in Ref. 2. The changes are 
for the most part refinements in analysis or program logic. A few minor 
errors in the programs have been found and corrected. The most extensive 
addition is the improvement of the use of tabular values for burning sur- 
face area during tailoff in the Monte Carlo program. This had been 
accomplished earlier for the design analysis program. 



II. TRF. MONTE CARLO EVALUATION OF PERFORMANCE PARAMETERS 


The ability to predict the thrust imbalance between a pair of SRMb 
firing in parallel, such as on the Space Shuttle or Titan IIIC, as a 
function of time is essential to the total vehicle system design. To do 
this with a reasonable deg .tie of accuracy without large scale testing 
programs has obvious economic advantage. As has been discussed in Refs. 

1-3, the Monte Carlo technique provides a reasonably accurate and economical 
means of predicting the performance variations of pairs of SRMs. In the 
present work, the utility of the Monte Carlo program has been extended to 
permit further evaluation of performance parameters. The Monte Carlo tech- 
nique is used to generate, as in Refs. 1-3, the performance of a theoretical 
population of motor pairs. Provision has now been made, however, for the 
data generated by the Monte Carlo program to be in the form of a magnetic 
tape containing the thrust imbalance versus time data for each SRM pair. 

The changes to the Monte Carlo program necessary to produce the tape are 
listed in Section VI of this report. 

The tape generated by the Monte Carlo program may be analyzed with 
respect to certain critical performance characteristics by the use of two 
new computer programs described in this section of the report. Such an 
analysis is extended for a sample case to a separate investigation of 
statistical correlations for a number of sets of pexfoimance parameters. 

As is discussed, one particular correlation may have important general 
practical application. 

The Imbalance Analysis Computer Programs 

Because of the various possible applications for the statistical 
analysis of the thrust imbalance versus time data and to maintain a certain 
amount of generality, the analysis was done in two parts. 

The objective of the first part of the analysis w.-’c, given any time 
interval during which the motors are firing, to determine the mean value and 
the standard deviation about this mean of the absolute value of the maximum 
thrust imbalance and of the time within the prescribed interval when this 
thrust imbalance occurs. This is calculated using a data search routine 
which is compatible with tl»' data tape generated ty the Monte Carlo program. 
The same calculations were also made for the time rate of change of thrust 
imbalance within the prescribed time interval. 

The computer program for the analysis of imbalance data during specified 
time intervals is listed in Appendix A. A sample of the computer output is 
shown in Table II-l. In addition, computer generated plots may be made by 
superimposing curves representing the thrust imbalances and thrust imbalance 
rates for each motor pair of the population. Tills feature is quite useful 
in deterndning the limits and general behavior of the imbalances for a 
particular population of motor pairs. An example of this graphical output 
is given in Figure IT-1. The latter example is based upon a tape of the 
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Table II-l. Sample printout for tine interval imbalance analysis computer program 


this is time interval number i 

THERE ARE 3d SETS OR DATA FOR THIS TIME INTERVAL 

THIS TIME INTERVAL BEGINS AT G.O SECS ANO ENOS AT 100.00 SECS 


THE ABSOLUTE VALUE OF 

this imbalance occurs 

the 

AT 

MAXIMUM THRUST 
7. 62 SECS 

IMBALANCE 

FOR 

MOTOR 

FAIR 

NUMBER 

1 

IS 

l.SSObl 

Ob 

LBF 

the ABSOLUTE VALUt OF 
this imbalance occurs 

THE 

AT 

maximum thrust 

B « VV SECS 

IM8AL ANCE 

FUR 

MOTOR 

FAIR 

NUMBER 

2 

IS 

2.A190E 

OS 

LBF 


I 

U 

the absolute value Or Th£ MAaInuN ThRUST IMBALANCE FOR MOTOR FAIR NUMBER SO IS 7»B>bOE OS LBF I 

THIS IMBALANCE r CCU MS AT 1.97 SECS 

THE ABSOLUTE VALUE OF THE MAXIMUM THRUST IMBALANCE DURING THIS TIME INTERVAL IS S.bbBBI 04 LBF 
THIS IMBALANCE OCCURS AT 8.0* SECS 

TMl MEAN VALUE UF 

the absolute VALUE OF THE maximum ThRUST IMBALANCE DURING THIS TIME INTERVAL IS 1.0S57E 0« LBF 
THE STANDARD OEVIATICN CF 

the ABSOLUTE VALUE OF The MAXIMUM THRUST IMBALANCE DURING THIS TIME INTERVAL IS 7.4TSSE OS LBF 
the MEAN VALUE OF the TIME FOR THIS IMBALANCE IS 22.89 SECS 


The STANDARD DEVIATION uF The TIME EUR THIS IMBALANCE IS SB. lb SECS 

THE ABSULUTE value UF the MAXIMUM THRUST IMBALANCE RATE FOR MOTOR FAIR NUMBER 1 IS 7.00bbC Ob LBF /SEC ' 

this imbalance rate cccurs at o.23 sics 

the ABSULUTE VALUE UF THE MAXIMUM THRUST IMBALANCE RATE FOR MOTOR FAIR NUMBER 2 IS 0.S7B0E 0) LBF/SEC 
This INjALANCC RATE OCCURS AT 0.22 SECS 

THE ABSOLUTE value of THE maximum Thrust IMBALANCE RATE FOR MOTOR FAIR NURBER S IS i.bObll Ob LBF/SEC 
THIS IMBALANCE RATE OCCURS AT 0.2b SECS 


THE ABSOLUTE VALUE OF the MAXIMUM thrust IMBALANCE RATE FOR MOTOR FAIR NUMBER 


b IS 


S.237SE Ob LBf/SCC 
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Figure II— 3 . 


CalConp plot of thrust imbalance versus time from 
0 to 100 sees, for 30 pairs of theoretical Titan 


IIIC SRHs . 
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performance of 30 pairs of Titan IIIC SRMs prepared using the input values 
described in Table 2 of Ref. 1. 

Any number of intervals may be specified for analysis. Also, similar 
data is obtained as output for the absolute value of the maximum impulse 
Imbalance and the times at which they occur. The absolute impulse imbalance 
is defined by the relation, 

t 

(4I t>Ab s'J |F 1- F 2 |dt HI - 1) 

o 

where Fj and F 2 are the thrusts of the two motors of a pair. The maximum 
absolute impulse imbalance clearly occurs at the last tabulated time point 
within the interval. 

The objective of the second part of the analysis was, given any time 
during which the motors are firing, to determine the statistical limits of 
the thrust imbalance and the first time derivative of the imbalance. In 
this analysis these limits are taken to be +K 2 S , where s is the second 
moment of the distribution about an assumed zero mean and ^ is the two- 
sided tolerance factor determined by the size of the sample being considered 
and the specified probability that a specified percentage of the distribution 
will be included within the limits. A separate data search routine is 
used to make these calculations and the results are presented in the form 
of computer generated plots. These plots determine a statistical imbalance 
envelope for the population. Sample plots are shown in Figures II-2 and 
II-3 based again on the 30 Titan IIIC pairs. Table II-2 illustrates the 
computer printout for this data. The computer program used for the time 
slice analysis is listed in Appendix B. 

The second program also generates data similar to that described above 
on impulse imbalance and absolute impulse imbalance. In the case of the 
absolute impulse imbalance the upper limit is taken as x + K 2 S where is 
the one-sided tolerance limit, s is the sample standard deviation, and 
x is the true mean because the + ^Sq limits would have little significance 
in this case. 

In order to establish the validity of this analysis, a comparison 
was made between a sample of 30 Titan IIIC SRM pairs obtained from the Monte 
Carlo program and data furnished by MSFC on 21 Titan IIIC flight tests. 

Figure II-4 shows that good agreement was obtained. Another comparison 
was made between 30 Space Shuttle type SRM pairs from the Monte Carlo 
program and the contractor's preliminary predicted imbalance envelope (Ref. 

7) during web action time. This comparison also showed good agreement with 
regard to the shape of the envelope as can be seen in Figure II-5. In- 
dividual SRM pair results are shown in Figure II-6 for the period 0 to 105 
seconds for the Monte Carlo evaluation. The contractor's prediction was 
based on the estimated 3 sigma spread between matched motors considered 
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Figure II-2. Tolerance limits for thrust imbalance from 0 to 100 
secs, for 30 pairs of theoretical Titan II1C SRJis 
(Probability 0.90 for 99.7% of the distribution). 
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Table II-2. Sample printout for time 
program. 


TMcRE Ml 30 SETS Of MOTOR PAIR OAT* 
the one S13EU A FACTOR (ATI FOR THIS SAMPLE SUE 
THE TOO StOEO K FAC I UR IK it fCH IMIS SAMPLE SUE 
The RESULTS ARE CALCULATED at 34 time slices 


THIS TIME SLICE MAS TAKE* AT Q.SO SECS 
the ♦ o* - R<2*stCMA limit about a zero mean for 
THE ♦ OR - R2*SICMA LIMIT ABOUT A ZERO MEAN FOR 


THIS TIME SLICE MAS TAKEN AT 5.30 SECS 

THE ♦ OR • *2*SICHA LIMIT ABOUT A EERO MEAN FOR 

THE * OR - K2*S!GMA LIMIT ABOUT A ZERO MEAN FOR 


THIS TIME SLICE HAS TAKEN AT 10.00 SECS 

THE * UR - KZ*S I CNA LIMIT AbOUl A ZERO MEAN FOR 

THE ♦ OR - K2*SIG*A LIMIT ABOUT A ZERO MEAN FOR 


this time slice mas taken at is. jo secs 

THf ♦ OR - KZ*S1CMA LIMIT ABOUT A ZERO MEAN FOR 
THE ♦ OR - RZ*S ICMA LIMIT ABOUT A ZERO MEAN FOR 


THIS TIME SLICE MAS TAKEN AT 20.00 SECS 
The ♦ OR - K2*SluMA limit about a zeru mean fur 
The » OR - K2*$!GMA LIMIT AoOUl A ZERO NtAN FOR 
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Figure II-4. Comparison of Monte Carlo tolerance limits tor 30 pairs of Titan III-C SKMs with 
flight test analysis for 21 pairs (K -3.68 and 3.86 respectively; probability 
0.90 for SS.7% of the distribution). 
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Figure II-5. Comparison of contractor's thrust imbalance during web 

action lime analysis (Ref. 7) with tolerance limits obtained 
by the Monte Carlo program for 30 space shuttle type SRM pairs 
(0.99 probability for 99.7% of distribution). 
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possible In throat erosion rate, delivered specific impulse, burn rate, 
and propellant weight. The combination of the variations produced a 2.0% 
burn time differential which is high with respect to the 0.85% differential 
predicted by the Monte Carlo program. We attribute the difference in the 
two results primarily to the much higher burn rete difference (C v *0.42%) 
used by the contractor with respect to that used in the Monte Carlo program 
(C v =0.06%). 

Correlation of T hrust Imbalance with Thrust Imbalance Rate 

For possible use in systems performance analysis, the correlation 
between thrust imbalance and its first time derivative was investigated at 
various operating times. This was done for both algebraic and absolute 
values of the variables using the data from the Monte Carlo program for the 
30 pairs of theoretical Titan IIIC SRMs. At each of the 20 time points 
considered the correlation coefficient r c between the variables was de- 
termined. The results are given in Figure II-7. Also, a linear regres- 
sion line was established for each time point. The slopes and intercepts 
of the linear regression lines are listed in Table II-3. 

Examination of che results shows that there is a better correlation of 
the algebraic values than the absolute values in the region from 10 to 80 
seconds, but the correlation coefficient for the algebraic values changes 
sign shortly before web time (107 seconds) . Both the algebraic and abso- 
lute values show poor correlation during the 15 seconds preceding and 
excellent correlation for the 15 seconds following web time. The varying 
characteristics of the regression lines at the various time points prevent 
further generalization of the results. 

It is possible that better correlation might be obtained by examining 
the relationship between thrust imbalance at one time and the thrust im- 
balance rate at a slightly earlier time. Such correlations have not been 
attempted because the potential applications of such works have not developed 
as originally anticipated. However, the analysis has suggested correlations 
between other variables which may have important consequences as discussed 
below. 

Correlation of Thrust Imbalance with Nominal Thrust Rate 

Investigation of the correlation between thrust imbalance and nominal 
thrust-time trace characteristics have revealed some interesting results. 

The mean of the absolute values of thrust imbalance at a given time was 
found to have a correlation coefficient of 0.976 with the absolute value of 
the nominal slope of the thrust-time trace at the same time. This is 
based on analysis of the previously discussed Titan IIIC Monte Carlo data. 

A similar correlation of the standard deviation of the absolute values of 
thrust imbalance with the nominal slope of the thrust-time trace gives a 
sample correlation "ocfficient of 0.954. Figures II-8 and 9 are scatter 
diagrams for the variables. When the thrust imbalance was correlated with 
both thrust slope and thrust, the correlation was not significantly im- 
proved. 
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Figure II- 7. 


Correlation coef f icionts between thrust imbalance 
and thrust imbalance rate at various times during 
operation for 30 pairs of Titan IIIC SRMs . 



-lb- 


Table II-3. Results of correlation and regression analysis of Monte Carlo 
thrust imbalance and thrust imbalance rate for Titan IIIC SRMs 
(AF - b AT + c) . 


Algebraic Values Absolute Values 


Time (sec) 

b 

c 

r c. 

b 

c 

*c 

0.5 

-0.6900 

-179 

-0.0774 

1.1291 

8114 

0.1590 

10. C 

-74.5114 

-2895 

-0.8109 

50.6067 

3871 

0.4905 

20.0 

-59.8190 

-1872 

-0.9242 

56.6595 

528 

0.7763 

30.0 

-7.9085 

-3215 

-0.8043 

3.3994 
. . . 

4329 

0.4001 

40.0 

-70.3361 

-241 

-0.9054 

6.0127 

1434 

0.7739 

50.0 

-90.7287 

-820 

-0.8476 

67.9270 

1960 

0.6126 

60.0 


-463 

-0.5535 

101.9460 

4138 

0.3736 

70.0 

^|ppPP| 

171 

-0.8339 

88.4369 

1603 

0.6678 

80.0 

-54.3596 

-412 

-0.8092 

44.2926 

1659 

0.6154 

90.0 

-46.6494 

-1049 

-0.5239 

46.2998 

2862 

0.5319 

100.0 

-53.3827 

-992 

-0.5424 

47.8907 

1028 

0.4708 

10 ' r 

-4.5988 

-232 

-0.2110 

3.6304 

554? 

0.2421 

ic : 

-0.0436 

1187 

-0.0907 

-0.0171 

4989 

-0.0526 

103.0 

0.5521 

3362 

0.7057 

0.5111 

4803 

0.6596 

104.0 

0.7497 

13709 

0.7070 

0.3391 

22439 

0.3124 

105.0 

1.2287 

3778 

0.4265 

0.8371 

42578 

0.4449 

106.0 

-27.8496 

2534 

-0.9731 

26.9581 

2929 

0.9324 

110.0 

-21.2082 

1279 

-0.9858 

21.1880 

273 

0.9651 

120.0 

-7.4078 

-141 

-0.9863 

7.4302 

-123 

0.9670 

125.0 

-0.5615 

3520 

-0.3174 

0.2900 

4899 

0.1642 
























































































































MEAN OF ABSOi TE VALUES OF THRUST IMBALANCE (lbf x 10~ 3 ) 



Figure II- 8 . Scatter diagram: mean of absolute values of thrust imbalance for 30 

theoretical Titan IIIC SRM pairs versus the absolute value of mean 
time rate of change of thrust for the 60 motors. 




STANDARD DEVIATION OF ABSOLUTE VALUES OF THRUST IMBALANCE (lbf x 10 



Figure II-9. 


Scatter diagram: standard deviation of thrust imbalance for 30 

theoretical Titan IIIC SRM pairs versus the absolute value of the 
mean time rate of change of thrust for the 60 motors. 
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The nominal slopes of the thrust-time trace were determined by 
computing the mean of the slopes for the 60 SRMs at each time point. 

It should be noted that the 20 time points on which the latter cor- 
relations were based were not randomly selected; they were selected to 
cover all portions of the trace for the previous analysis of the relation- 
ship between thrust imbalance and imbalance rate. Although the time points 
were not selected with a knowledge of the values of the present variables 
of interest, the manner of selection introduces a question on the random 
nature of the variables (thrust imbalance and thrust slope), so that further 
investigation of the correlation should be accomplished before attempting 
additional statistical analysis. However, it appears that there is a 
strong relationship between the variables in question. 

The relationship between the thrust imbalance and the slope of the 
thrust-time trace has at least two important potential applications. The 
first can be accomplished only if the regression analysis proves to be 
general, that is, it must be shown that it is reasonable to expect that the 
same regression line fits all SRMs or at least members of a family of SRMs. 
If it does, the thrust imbalance predictions could be based directly on 
the slope of the trace of a new SRM without Monte Carlo evaluation. Also, 
although it has been shown (for Titan IIIC's, at least) that the Monte 
Carlo program gives reasonable representation of the SRM pair imbalance, it 
would be a logical step to test the correlation directly against actual 
performance data. This suggests a second application - that of using the 
regression analysis to establish imbalance scaling relationships between 
different SRMs. The regression analysis for the Monte Carlo program would 
be compared to that for the actual moror to establish the scaling relation- 
ship. The result would then be applied to a regression analysis of a Monte 
Carlo evaluation for a new SRM (e.g., the Space Shuttle) to establish the 
imbalance limits for the new SRM pairs . This could be done whether or 
not the same regression line fics the two types of SRMs. Investigation of 
these applications is beyond the scope of the current research. 



III. EFFECT OF GRAIN DEFORMATION ON' INTERNAL BALLISTICS 


The circular-perforated (c. p . ) portion of the grain of an SRM of the 
Space Shuttle type may deform as much as 1.2 inch in the radial direction 
under pressurization from the chamber gases. At a 30-inch bore radius, 
this results in as much as a 42 increase in the initial burning perimeter 
of the grain bore. This can affect the chamber pressure and the apparent 
burning rate of the propellant. The effects will decrease as burning of 
the propellant web progresses because the displacement of the burning 
surface relative to its unpressurized position will decrease. The ability 
to quantitatively assess these particular effects of grain deformation on 
internal ballistics should improve predictability of performance of indi- 
vidual motors, especially for the first motor of a new design where the 
"scale factor" on burning rate is uncertain. 

Analysis of Grain Deformation Effects 


The underlying hypothesis on which the ballistic effects of the grain 
deformation is based is that at. fixed pressure the regression rata r of 
the pressurized propellant surface just beneath the heat-affected c 
zone is independent of the state of strain. To show that this is plausible, 
first consider the burning perimeter of a sector of solid propellant. The 
burning perimeter consists of a thin zone (liquid and/or solid) in which the 
physical properties are degraded to a state where the propellant develops 
negligible shear stresses. Therefore the regression rate of this zone (r ) 
should be independent of the precise state of strain calculated just beneath 
the degraded zone which will be somewhat thinner than the solid phase heat- 
affected zone. The concept needs experimental verification, but appears 
consistent with the known features of the solid-gas transition region. 


In the present analysis, it is assumed the thermoelastic coupling 
discussed in Section IV on this report is absent. Thermoelastic coupling 
may possibly alter the depth of and the temperature distribution within 
the heat-affected zone and thereby make burning rate sensitive to strain. 
However, as discussed in Section IV, it is expected to be significant only 
during highly transient conditions; i.e., during ignition and very rapid 
depressur izat ion . 


If we now consider a control surface (See Fig. III-l) separating the 
region of degraded propellant from the remaining propellant, a quasi-steady 
state mass balance across the heat-affected zone, which for the purpose of 
this analysis is considered of zero thickness, yields 


r 

c 


= r p /p 
mm pc 


(III-l) 


where p is the density of the degraded propellant and p is the density of 
the unheated propellant just outside the heat-affected zBne. It is apparent 
that r and p are pressure dependent but should not depend on the state of 
strain underneath the degraded surface. Also, when as in the principal case 
of interest here, Poisson's ratio v is close to 0.5, it is intuitively clear 
and may easily be shown that the density of the unheated propellant can be 
expressed as a function of pressure (and modulus) only since 


(l+e r ) (l+e 0 ) (l+c z ) = U-P(l-2v)/E] 3 


(III-2) 
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Figure. I1I-1. Sector of cross-section of solid propellant with 
and without pressurization and combust ion. 



Figure 1II-2. Grain geometry notation for analysis of solid- 

propellant strains. The insulation and liner may 
be considered a part of the propellant. 
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where e is strain, P is pressure, E is modulus o£ elasticity and 8, r and s 
refer to the tangential, radial and axial directions, respectively. Thus 
Eq. III-l should govern the burning rate of the propellant regardless of 
the state of strain and r is a function of pressure only (for fixed 
Initial grain temperature c and erosive burning characteristics) . 

Next, the effect of the deformation on internal ballistics is examined. 
We consider only a circular-perforated grain for which e fi is independent of 
the tangential coordinate at one axial position. The mass generated per unit 
time per unit length at the control surface separating the heat-affected 
cone of the propellant from the unheated propellant is given by 

m G /l - S o (l+e o ) Ppo r c /{ [l-P(l-2v)/E] [14u(T Gr -T Ref ) 1 > 3 (1II-3) 

where S is the burning perimeter, a is the linear coefficient of thermal 
expansion and the subscript o refers to the impress* rized and unheated 
(reference grain temperature) state of the propellant. The approximation 
of Eq. 111-2 has been used in Eq. III-3 to determine the effect of pressur- 
ization on the unpressurized density, and that density has been further 
modified for the small effect of thermal expansion between the standard 
temperature T , at which density is measured and the actual grain tempe- 
rature T . Ref 
Gr 

A second modification of customary ballistic analysis is the way the 
time At required to burn an increment Ay normal to the surface is calculated. 
The usual relationship is 


At = Ay/r (III-4) 

o 

where r is the burning rate deduced from strand burners or small ballistic 
test motors modified with a "scale factor" reflecting the estimate of the 
change in burning rate between the strands or test motors and the SRM under 
consideration. The comparable equation to Eq. III-3 used with Eq. III-4 is 

m _/£ = S p r (III-5) 

G p po o 

For the present we disregard any scale factors and devise an appropriate re- 
lationship for At in terms of Ay and r applicable to the deformed grain 
problem. First, note that the r shouSd be the rate at which the unpres- 
surized and unheated propellant regresses because the Ay is to be used for 
further computations of surface area as an increment of undeformed propellant. 
Next, it is assumed that the changes in length (Z) of the grain produced by 
pressure and thermal loading are negligible. As shown in Ref. 4, longitud- 
inal (axial) grain strain due to pressurization is on the average over the 
length of the grain approximately a half order of magnitude less than 
the corresponding tangential strains. The axial strains are negative while 
the tangential strains are positive. Although a quantitative assessment 
has not been made, when heating of the surface of the propellant is considered, 
the axial strain at the point of interest; i.e., just beneath the heat-affected 
zone, will become less negative because of the restraint of the unheated 
propellant to the expansion of the heat-affected propellant. ”'iis will de- 
crease the extent of any change in length of the propellant relative to the 
increase in perimeter over that due to pressurization alone. 
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Ref erring to Fig. III-l, ve obtain from a mass balance between deformed 
and undeformed propellant 


p S Ir ■ 
po o o 


P n , s l(l+e )r 
pc O 0 C 


(II1-6) 


or 


r Q - r c (l+c 0 )/{U-P(l-2v)/E][l+a(T Gr -T ref )]> 3 (III-7) 

The same result may also be obtained by direct combination of Eqs. I1I-3 and 
III-5. With this r Q , Eq. III-4 then applies to the deformed grain situation. 

Data to give r are obtained from ballistic test motors or strand 
burners. As previously noted, no scale factor is considered. Indeed, it 
is suggested that the grain deformation may account for a significant 
portion of the scale factor. Ihe strands are under essentially a state of 
hydrcstatic compression. The same is very nearly true for the typical 
ballistic test motor because of the rigidity of the relatively thick cases 
and the relatively thin webs. It should be recalled that r should be inde- 
pendent of the state of strain. However, the burning rate of small test 
motors (or strands) is deduced by dividing the undistorted web thickness 
(or strand length) by time. Therefore, to determine r , the regression rate 
of the propellant surface just beneath the heat-affectld zone of the SRM, 
the test motor (or strand) rate must be reduced by the multiplying factor 
(l-P(l-2v)/E], The r , of course, is subject to the usual laws governing 
its relationship with c initial grain temperature, pressure and velocity (erosive 
burning) . 

Simplified Determination of Grain Deformation 

For a first analysis of the effect of grain deformation on internal 
ballistics, a simplified model of the deformation was developed. We are 
concerned for the present only with the modification of the internal burning 
surface of a circular-perforated grain for which it has been shown that the 
important parameter is the tangential bore strain e 0 . The will change 
as burning progresses and also will vary in general with axial position. 

Smith (Ref. 4) working under the direction of this project has shown that 
as an approximation the may be considered independent of axial position. 
Smith extended the method used by Vandenkerckhove (Ref . 8) for calculation of 
grain stresses in case-bonded propellant with rigid motor cases to the 
calculation of propellant strains in non-rigid motor cases. His results for 
tangential strains at the bore are given by 

E e n /P. = {[1+m 2 + v (1-m 2 - 2m 2 v )]/(l-m 2 )} 
p 0 i p p 

+ (P /P,M [2(v 2 -l)/(l-m 2 )] + E R v v /E t } (III-8) 

ei p pepccc 

where P £ is the radial stress at the propellant-case interface as given by 

P /P. = [2m 2 (l-v 2 ) / (1-m 2 ) ] / { [E R (1-v v )/E r ) - v (1+v ) 
ex p pecpcc PP 

+ (1-Vp 2 ) (1+m 2 ) / (1-m 2 ) } 


(II1-9) 
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Here the geometric notation is given by Fig, III-2, is internal pressure, 
v is Poisson's ratio and E is modulus of elasticity. The subscripts p and 
c refer to the propellant and motor case, respectively. Details of the 
analysis are given in Ref. 4. 

Smith compares results for tangential, radial and axial strain with a 
more exact solution obtained using a computer program devised by Brisbane 
(Ref. 9). The two solutions for the tangential strain are illustrated in 
Figs. Ill- 3a and b for a 140-inch diameter motor segment 300 inches long for 
various web thicknesses. The results are given in non-dimensional form and 
are applicable to motors of various sizes with the same ratios T w /R e , r c /R e 
and 1/Rgand material properties. Study of Figs. III-3a and b shows that the 
simplified solution, which assumes the same strain at all axial positions, 
in general, overestimates the strain with respect to the more exact solution. 

The worse case is at the largest web thickness to outside radius ratio (TAU/RO) . 
Here, the strains calculated by the two methods are within 30Z over two-thirds 
of the grain length and the comparison improves rapidly with decreasing web 
thickness (increasing distance burned) . It should be noted that the strains 
calculated by both Refs. 4 and 9 are those due to pressure loadings only. 
Expansion of the propellant in the heat-affected zone will tend to produce 
additional tension of the solid surface just outside the heat -affected zone 
(the surface of interest) . This will compensate somewhat for the apparently 
excessive strains calculated by Eq. III-8. 

We feel that the approximate analysis is adequate for a first approach 
to analysis of the described effect of deformation on internal ballistics. 

It will permit rapid calculation of the bore strains for various configura- 
tions, material properties and distances burned. Additional investigation 
performed by the authors under this project using the program of Brisbane 
shows that tangential bore strain, the quantity of primary importance in the 
analysis, is changed by less than 5Z by solid-phase thermal gradients so 
that thermal effects may be neglected in this first analysis. The reason 
the thermal gradient changes the tangential bore strain so little is that 
the gradients are confined to only a very thin zone at the surface. The 
bulk of the strain is due to pressurization of the entire propellant web. 

Modification of the Design Analysis Computer Program 

The design analysis computer program of Ref. 2 has been modified to 
permit computation of the internal ballistics with the deformed grain. The 
program is listed in Appendix C along with instructions on preparation of 
new input data cards and a sample problem solution. 

Bore deformation of c.p. grains or grain segments is evaluated by the 
simplified approach discussed earlier. Only deformation of the internal 
burning surface of c.p. grains or segments is considered. The calculation 
of the deformations of the ends of c.p. grains and star grains is a formid- 
able task. Any solution would appear to be inconsistent with the simplified 
analysis of deformation and internal ballistics being presented and well 
beyond the scope of the present research. The design analysis treats the 
ends of c.p. grains and star grains with an unmodified burning ratio. This 
is accomplished by computing a separate burned increment Ay' applicable to 
c.p. grain ends and the entire star grain. The Ay* (computer symbol YETA) 
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Figure III-3a. Comparison of approximate solution for tangential 
bore strain (Kef. 4) with finite element solution 
method of Ref. 9. 
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is related to the Ay applicable to the c.p. grain lateral surface by the 
equation 


Ay* - Ay{[l+P(l-2v p )/E p ][l+a(T Gr -T ref )]} 3 /(l+E e ) CHI-10) 

The Ay' thus calculated gives the regression of the star grain surfaces and 
ends of the c.p. grain in the same time increment that the c.p. lateral 
surface regresses a distance Ay. The modification of the bore surface of 
the c.p. grain due to deformation is accounted for in the mass generated 
equation by the use of the program variable XETH which is the fraction of 
the total burning bore surface associated with the c.p. grain. The factor 
(1 + XETH*ETHETA) , where ETHETA ■ Gq, multiplies the undeformed total burn- 
ing bore surface area when the deformation option is elected by setting 
ISOl, The program will calculate the internal ballistics without grain 
deformation when IS0=0, The program changes that have been made in the 
design analysis of Ref. 2 to permit the alternative computations are identi- 
fied by the symbol > in the left-hand margin of the listing. Also, the 
modifications of the design analysis program discussed in Section VI of 
this report have been incorporated into the program listed in Appendix C. 

In modifying the program, it was found convenient for the purposes of 
minimizing the number of modifications and for preserving the ability of 
the program to make calculations with no deformation to change the usage 
of SUMAB in the program. Before the modification (or presently when IS0=0) 
SUMAB represents the surface area at a given distance y burned normal to 
the entire surface. With the modification (when IS0=1) because of the 
difference between y and y* the SUMAB represents the area at a distance y 
burned over a part of the surface and y f over another part. This creates an 
error in the calculation of the weight of propellant burned (WP2) which is 
obtained by integration of SUMAB with respect to y alone. Therefore the 
calculation of WP2 should be disregarded when IS0=1. To calculate WP, the 
arithmetic average of WP2 and the weight of propellant burned (WP1) calcu- 
lated from the mass discharged rate, the program now sets WP2=WP1 when ISO=l. 
Comparisons of weights calculated by the two methods can still be made by 
comparing WP1 calculated with IS0=1 with WP2 calculated with IS0=0. The 
latter is still a valid calculation for total weight of propellant burned 
even with grain deformation. Also comparisons have been made for the two 
sample cases discussed next and the WP1 is not significantly different 
for IS0=0 and IS0=1. 

Comparison of Theoretical and Experimental Results 

To test the validity of the hypothesis regarding the effect of grain 
deformation on internal ballistics, the modified program was used to pre- 
dict the performance of two entirely different rocket motors and the 
results compared with actual performance data. In each case the burning 
rate used to predict the performance was that obtained from ballistic test 
motors or strand Durners; no scale factor was applied. 

The first comparison made was for a Titan IIIC/D configuration. The 
basic input parameters used for the prediction are those given by the mean 
values in Table 2 of Ref. 1 with the following exceptions: The propellant 

density had been changed to 0.06325 lbm/in 3 because additional data on this 
has become available. The strand burning rate coefficient obtained from a 
number of batches of propellant loaded into the same 5RM (No. D5) for which 
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the comparison is to be made was used for the prediction. The burning rate 
coefficient used (0.06466 with no scale factor) was deduced from Ref. 10. 

Five burning slot faces were assumed in the program calculations. A table 
of tabular values was prepared so that the burning surface areas calculated 
by the unmodified program more nearly match the known geometric surface 
versus distance burned characteristics of the Titan II1C/D configuration 
over web action time. Aft-end closure propellant surface representation 
was included in the tabular values. 

Propellant and motor case properties used in the analysis of the de- 
formed propellant grain are given on Fig. III-4. The figure is a comparison 
of the program results for aft-end stagnation chamber pressure versus time 
with and without deformation and with the results of analysis of actual test 
data given for Titan IIIC/D motor number D5, in Ref. 10. This particular 
motor was selected because its performance is more or less typical of the 
Titan IIIC/D and because of availability of strand burning results on the 
same propellant loaded into the large motor. 

A similar comparison of theoretical and experimental results was made 
for the Castor TX354-5 SRM. This is a 31-in. dia. motor with an entirely 
c.p. slotted grain which is described in detail in Ref. 11. The SRM is 
also described briefly in Ref. 12 where input parameters applicable to an 
early version of the design analysis program are identified. A complete 
listing of input values for the modified design analysis is given in Table 
C— 1 of the present report where the Castor motor is used as a sample problem 
to illustrate use of the program. The input values differ from those used 
in Ref. 12 because of the program changes that have been made and because 
we have improved our representation of the grain geometry, particularly with 
respect to use of tabular values. Again, tV. burning rate coefficient used 
has no scale factor between the Castor motor and the small ballistic motor 
data from which the coefficient was deduced. Vacuum thrust results are 
compared (See Fig. III-5) for the Castor motor. The vacuum thrust results 
for the actual motor were based on analysis of static test results and are 
tabulated in Ref. 11 as typical performance values. The strand burning 
rate data is also typical but are not otherwise identifiable with the 
specific large motor. 

Analysis of the two comparisons shows that the theoretical predictions 
are substantially improved by consideration of the deformations. Additional 
comparisons are needed to confirm the efficacy of the new model. Modifica- 
tion of other parameters can also make improvements in prediction possible. 

A prime example is the use of the burning rate scale factor. However, the 
scale factor is most often applied without specific theoretical justification 
and usually is only accurate after data on actual SRM firings has been ob- 
tained. On the other hand, the deformation effect has a basis in theory. If 
it can be further validated, it will not only by itself improve prediction 
capabilities but also clear the way for more accurate assessment of other 
specific ballistic phenomena such as erosive burning and nozzle throat erosion. 
For example, consider the discrepancy between the prediction with grain defor- 
mation and the test results for the Titan IIIC/D. It is obvious that the 
areas under the two traces are not the same. This is suggestive of a differ- 
ence between the nozzle throat erosion rate used in the theoretical prediction 
(time independent, pressure and size dependent) and the actual erosion rate. 

The application of a scale factor to the theoretical burning rate coefficients 



AFT-END STAGNATION CHAMBER PRESSURE (psi) 



20 40 60 80 100 120 140 

TIME (Seconds) 


Figure III-4. Comparison of theoretical and experimental results for a Titan IIIC/D SRM 
with no scale factor on strand burning rate. 





Figure III-5. Comparison of theoretical and experimental results for a Castor TX 354-5 SRM with 
no scale factor on small ballistic test motor burning rate. 
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vould not improve the area discrepancies and would mask the true throat 
erosion and/or buildup phenomena. 

Variations in the propellant properties used in the deformed grain 
program can affect the results. Likewise variations between SRHs of a 
pair firing in parallel can contribute to thrust imbalance. As shown 
in Ref. 4, the effect of propellant modulus is especially important. 
However, additional confirmation of the deformed grain model, preferably 
Including direct experimental verification, should be obtained before 
the effect is Incorporated into the Monte Carlo program. 

Equations III-3 and III-7 which give the basic burning rate theory of 
the model may also be coupled to more rigorous descriptions of the grain 
deformation and/or more rigorous models of the other aspects of internal 
ballistics than those used in the modified design analysis program. 



IV. THERMOELASTIC ANALYSIS 


The problem to be discussed In this section is the effect on the 
solid propellant burning rate produced when heat is generated by mechanical 
deformations of the solid propellant. The basic theory and the numerical 
solution procedure used in the thermoelastic analysis as well as the 
results which have been obtained are presented. 

The basic formulations and concepts associated with the coupling 
between thermal and mechanical loadings can be traced to the work of 
Duhamel (Ref. 13). Other early investigators include Lord Kelvin (Ref. 14) 
and Joule (Ref. 15) . It was Duhamel who developed the relationships used 
to express the distribution of strain in a body subjected to differential 
heating. 

Basic Theory 

Consider a small element of an isotropic elastic solid subjected to 
an arbitrary stress which is removed from its surroundings and subjected 
to a temperature change, AT. The additional strain in the element is given 
by the tensor BAT6j.,/3, where 8 is the coefficient of volumetric expansion 
and is the Kronecker delta. From this it follows that if the strain 
tensor is originally given by ejj , then the new strain tensor would be 
e^j - f?AT(5ij/3 and this would be the strain used in the constitutive rela- 
tion. The effect of a nonuniform temperature distribution becomes equiva- 
lent to that of a body force per unit volume, given by - (B/n)V(AT), n 
being the compressibility and V the gradient operator. Thus, once AT is 
known from the solution of the heat conduction problem, the system of 
equations for the displacement field is fully defined. As can be seen, 
this approach postulates the effect of the temperature gradients on the 
state of strain, but assumes that the heat transfer is totally unaffected 
by the state of strain. This result is upheld rigorously when the medium 
is in mechanical and thermal equilibrium, but its application to time 
dependent problems is not satisfactory. However, Duhamel and Neumann do 
suggest that a tern proportional to Sr^/31, the rate of change of dilatation, 
be included in the heat conduction analysis. The present work investigates 
this term's effect on a thermoelastic analysis, and in turn, on the pro- 
pellant burning rate. 

Much effort has been expended in recent years to obtain solutions for 
problems involving thermomechanical coupling. This has been due, at least 
in part, to the necessity for structural analysis of bodies subjected to a 
nonisothermal, transient tenperature distribution. One of the better ex- 
amples of a body subjected to this type of loading is the solid-propellant 
rocket motor. The general analytical formulation of the problem is well 
known and can be found in various publications, of which Ref. 16 is typical. 
The problem to date has been obtaining solutions for problems whose geometry 
and loading conditions, both thermal and mechanical, are sufficiently general 
to be of practical interest. 
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Several approaches appear to be suitable for analyzing problems in- 
volving thermomechanical coupling. The first possibility is that of ob- 
taining analytical solutions of the governing equations. There have been 
several such solutions in the literature in recent years. However, the 
approach to now has been severely limited as to choice of geometry, loading 
and boundary conditions. The second possibility utilizes the finite element 
method and entails the approximation of the full set of the governing, 
coupled, nonlinear partial differential equations by a set of nonlinear 
algebraic equations which are then solved numerically. This finite element 
method is considerably more flexible than the first approach with regard to 
the variety of geometries, loading and boundary conditions which can be 
treated and has been applied to several classes of problems. However, this 
nonlinear finite element method also appears to be somewhat limited due to 
the complexity of solving large numbers of simultaneous nonlinear algebraic 
equations. Also added complexity is introduced if it is desired to investi- 
gate the effect of nonlinear material properties (except for materials 
which obey relatively sisple mathematical laws), arbitrary loadings, and 
boundary conditions. 

A third possibility is to approximate the transient nonlinear solution 
by & series of linear solutions. In this approach the variable, time, is 
included in what basically can be described as a time stepping procedure. 

At each time step, the initial thermal and mechanical states are established 
based on the conditions which are specified either by a direct input or by 
the conditions which existed at the end of the preceding time step. The 
solution then proceeds by obtaining the temperature distribution in the body 
for the given thermal loading. The thermal load may include applied tem- 
perature distributions on the external boundaries and/or at points within 
the body, heat flux through the boundaries, and convective heat sources (or 
sinks) distributed throughout the body. The internal heat sources are 
particularly important since they are used to account for the coupling 
phenomena. The solution then proceeds by using the computed temperature 
distribution as an input to the stress analysis portion of the computational 
procedure. Determination of the deformations, strains and stresses due to 
the combined loading is then made. Once this is done for a particular time 
interval, the coupling term, which is a function of the material properties 
(which may be functions of temperature), the local temperature, and the local 
strain rate, can be evaluated. With this information the computation then 
proceeds to the next time step where the coupling term is included as an 
internal heat source in the heat transfer analysis. 

This third approach, the one used in this research, appears to eliminate 
most of the problems associated with arbitrary time dependent loadings, 
temperature dependent material properties, and also some of the numerical 
computational problems associated with the nonlinear finite element method 
described above. This reduces the numerical computation difficulty because 
the problem is now formulated as a set of linear algebraic equations. The 
main objection to this approach is obviously the linear modeling of a non- 
linear problem. However, if the mathematical modeling of the geometry. 



-34- 


the loading, the material properties, and the time step is done carefully, 
this method is sufficiently reliable for the solution of a considerably 
larger group of problems than was heretofore practical. 

The work presented here is restricted to axisyimnetric bodies. This 
does not imply a restriction of the method, but is made because of the 
applicability to the problem at hand. 

Problems in thermoelasticity require not only the solution of the 
equations of elasticity, but in addition those equations necessary to 
describe the thermal state of the body being analyzed. In many problems 
the thermal state can be sufficiently defined by determining the distribution 
of temperature throughout the body and also any changes in this distribution 
with respect to time. Quite often as in the present work it is convenient 
to obtain the temperature distribution from the conservation of energy 
equation. The discussion which follows will be primarily concerned with 
the energy equation since it is the equation which will be used to determine 
the Influence of the thermoelastic coupling on the propellant burning rate. 

The elasticity equations are not presented because they are used in a more 
or less classical manner. 

As in Ref. 2, the conservation of energy equation for a control volume 
moving with a velocity r(t) with respect to an observer situated on the 
boundary (x=0) of the control volume and modified to account for volumetric 
heat release or absorption within the solid phase is 

3T/3t - r(t) 3T/3* + (X/pc) V 2 T - ToE e/[pc(l-2v)J (IV-1) 

where x is measured radially outward from the surface, e represents the sum 
of three translational strain components at the point under consideration, 
and the entire last term represents the heat dissipated by elastic deforma- 
tions . 

The two terms which are >f particular interest in the present analysis 
are the first and last terms on the right-hand side of Eq. IV-1. As mentioned 
above, the last term represents the contribution of heat from the elastic 
deformation of the body and can be calculated directly if the instantaneous 
local temperature is known as well as the time rate of change of the trans- 
lational strain components. Once the initial temperature and strains are 
prescribed, the temperature at any point within the control volume at any 
time is obtained from the solution of Eq. (IV-1). The strain rates can then 
be calculated at any instant of time once the mechanical load is prescribed 
as a function of time. Hence, with regard to its calculation, the thermo- 
elastic coupling term offers very little difficulty. 

On the other hand the first term on the right-hand side of Eq. IV-1 
requires that the instantaneous value of burning rate be calculated. One 
might initially propose that the burning rate r(t) could be easily obtained 
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from the pressure ?(t) which is required for the thermoelastic term by ex- 
tending the “standard" burning rate lav for steady-state burning to the 
transient problem as 


r(t) - a P(t)“ (IV-2) 


However, as pointed out in Ref. 5, the instantaneous burning rate r(t) may 
differ greatly from that obtained from a steady-state analysis due to a 
highly complex heat feedback mechanism between the various regions in and 
near the combustion zone. Therefore, it is necessary to construct a model 
to account for this thermal feedback mechanism in order to properly 
evaluate r(t). 

One such model for determining the transient burning rate developed by 
Zeldovich and Novozhilov, referred to as the Z-N model, is described in Ref. 
5. The present work uses this model to evaluate r(t) which is then in- 
corporated into the conservation of energy equation. As in Ref. 5, the Z-N 
method uses measured steady-state burning rates as a function of pressure 
and ambient temperature. This data is then used to construct the heat feed- 
back function from the gas to the solid in the proper instantaneous form to 
be used to study transient burning problems. The assumption which is made 
to allow use of steady state data is that the rate processes in the gas phase 
and at the propellant surface can be considered quasi-steady in the sense 
that their characteristic times are short compared to the pressure transient. 
The validity of this assumption is demonstrated in Ref. 5. What develops 
from this is that the steady state heat flux and the transient heat flux 
have the same functional form, which is written for the transient problem 
as 


X(3T/3x) , = r pc(T - T ) (IV-3) 

surface so 

where T s is the instantaneous solid surface temperature and T 0 is a 
fictitious temperature which exists infinitely far from the solid surface. 
Now, assuming the pyrolysis law to be of the Arrhenius type 


T - - E /[RinCr/A ) ] 
s s s 


(IV-4) 


and the empirical expression for r in terms of p and T Q to be given by 

r = a p n exp[a (T - T ,)) 
v r p o ref 1 


CIV-5) 
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vhere Op is the temperature sensitivity coefficient at constant pressure, 
which for this work is expressed as a linear function of pressure by 


o p - 1.1 x 10- 6 P + 6.9 x 10- 4 (IV-6) 

One can solve simultaneously Eqs. IV— 1 , 4 f and 5 to obtain the unknowns 

r, T and T . 
s o 

Numerical Solution Procedure 


The finite element method as used in this study is based on the work 
of Becker and Parr (Ref. 17). The basic program used for the finite 
element analysis is that of Brisbane (Ref. 9) and is primarily an extension 
of that given in Ref. 17. The method used is also quite similar to that 
given by Uilson and Nickell (Ref. 18). 

The use of the finite element method involves the division of the region 
of interest into a number of subregions which normally have a relatively 
simple geometric shape. Usually a polygon is used for two-dimensional prob- 
lems and quite often triangles are used. In this work the problem is axisym- 
metric and a circular ring element with a triangular cross-section is used. 

A typical element is shown in Fig. IV-1. The subregions, or elements, are 
connected at discrete points called nodal points. These nodal points are 
usually at the comers of the polygon as will be the case in the present work. 

For the heat conduction problem a simple relationship for temperature 
and its first derivative as a function of position within the element is 
assumed. One of the more common assumptions and the one made here is that 
both temperature and its first derivative are linear functions of position 
within the element. Hence, for the axisymmetric problem the temperature and 
its gradient can be written as 


and 


T(x,z) = + C 2 x + 0^2 (IV-7) 

9T(x,z )/3t = C 4 + C 5 x + C 6 z (IV-8) 


It should be noted that T(x,z) and 3T(x,z)/3t are considered to be evaluated 
at a particular time; e.g., t = t*, and hence, no explicit time dependence 
is expressed in the above equations. 

Once this approximation has been made, it is substituted into the basic 
field equations which describe the heat transfer process. The result is 
a set of linear algebraic equations, the unknowns of which are the temper- 
atures which exist at the nodal points. 



- 37 - 





-38- 


Basically the same procedure is used for the elasticity problem except 
that the displacements are approximated instead of the temperature. Another 
linear approximation is used for the displacements. This approximation is 
then substituted into the field equations for the elasticity problem and a 
set of algebraic equations is obtained for the nodal point displacements. 
These nodal displacements are used to evaluate the constants in the linear 
displacement assumption which is then appropriately differentiated to obtain 
the strains at a particular time. 

It is now appropriate to examine in more detail the conqmtational pro- 
cedure which is used and the manner in which the coupling between the thermal 
and mechanical loading as well as the burning rate term are introduced into 
the numerical computational process. The solution technique consists of 
solving separately the energy equation and the elasticity equations at each 
time ctep. Outputs from each are then used as inputs to the other, thereby 
forming a closed loop computational scheme. The calculation of the tem- 
perature dependent term which appears in the stress-strain relations is 
common to both the uncoupled and coupled solutions. It is calculated from 
the apparent element temperature which is obtained from a geometrical 
weighting of the four nodal point temperatures which are computed in the 
heat conduction part of the analysis. 

The computation of the thermoelastic coupling term which appears in 
the energy equation, Eq. IV-1, requires some knowledge of the deformation 
history of the body in order to compute the strain rate. For the work pre- 
sented here it is only necessary that the state of strain at time, t - At, 
and at time, t, be known, where At is the time increment between steps. 

This requires that the strain components for each element be retained at 
the end of each time step. The old values are then replaced by the new 
state of strain which is then used in the next time step. The initial 
strain has been taken to be zero, but no significant problems arise if a 
non-zero initial strain condition is used. For an isotropic body, the 
thermoelastic coupling term in Eq. IV-1 involves only the dilatational 
components of the strain tensor. The axisymmetric solution requires that 
three strain components be known for each element. However, the coupling 
term in Eq. IV-1 contains the sum of the three dilatational components so 
that only the trace of the strain tensor need be retained for each element. 
This results in a considerable saving of computer storage. The strain rate 
for each element can then be computed by taking the difference between the 
trace of the element strain tensor at time, t - At , and at timt., t, and 
dividing by the increment At. The resulting linear approximation is then 
transferred to the heat conduction analysis. The remaining variables in 
the coupling term involve the temperature and material properties of the 
element. The temperature is known from the calculations described above 
and the material properties are either prescribed constants or are given 
in tabular form as a function of temperature. If the material properties 
are given as constants, a direct substitution is made and the coupling term 
is evaluated directly. If the material properties are given in tabular form, 
a linear interpolation routine is used to obtain values not given in the 
table. 
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The term containing the burning rate requires that the Z-N model be 
used to evaluate the instantaneous burning rate r(t) based on the existing 
surface temperature at the end of each time step. The burning rate term 
is then added to the heat balance for each element. The heat flux at the 
surface is computed using Eq. IV-3 and this becomes the heat flux boundary 
condition for the next time step. 

The numerical example presented here is a simulated ignition of a 
large SBM. The motor is modeled as a straight circular perforated grain 
with an initial bore diameter of 60 inches and a web thickness of AO inches. 
The grain is encased in a 0.5-inch thick steel case. The motor is subjected 
to a time dependent Increase in pressure on the bore and end surfaces, 
but burning is allowed in a direction normal to the bore surface only. The 
variables in the analysis are assumed to be distributed symmetrically about 
the motor centerline. 

Finite element model 


It is well known that solid propellants are extremely good insulators 
and that the heat-affected zone of a burning propellant is of the order of 
about 10 mils. The finite element model is constructed to analyze this 
region with reasonable accuracy. A schematic of the finite element model 
is shown in Fig. IV-2, Finite elements with a radial thickness of 1 mil 
are used for the first 20 mils of the propellant web (Steady state tem- 
perature profiles show that this will more than adequately cover the heat- 
affected zone under most normal operating conditions). The remaining web 
is modeled using finite elements which have a radial thickness of approx- 
imately 10 inches. The steel case is modeled by two elements 0.25 inches 

thick. The vertical dimension of each element is 10 inches. The coarse- 
ness of the grid in the vertical direction is justified because the vari- 
ables under consideration (i.e., temperature, displacements, strains, etc.) 
do not have gradients in the vertical direction. To further refine the 
numerical calculations, each rectangular cross section element described 
above is subdivided into four triangular elements of the type illustrated 
in Figure IV-1. 

As the propellant regresses, a new finite element grid is constructed 
at each time step. The f irst 20 mils of the propellant is always modeled 

by elements of the size described above. The remaining elements, except 

those associated with the steel case, become decreasingly small as the 
web thickness diminishes. This is done to assure that the same degree of 
accuracy is maintained with regard to the propellant xn or near the heat- 
affected zone. 

Initial conditions 


The initial conditions are prescribed as: 1) no initial strains exist; 

2) an initial uniform pressure exists over the entire propellant surface; 

3) the initial temperature distribution corresponds to that which would 
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exlst if the propellant were burning under steady state conditions at the 
existing pressure. The equation for this distribution is 


• T(x) - [T s ( Pi ,T o )-T o ]ex P [- (xpcr/l)] + T q (IV-9) 

where T is taken to be the Initial propellant bulk temperature; 4) the 
initialneat flux at the burning surface is obtained using steady state 
values in Eq. IV-9. 

Boundary conditions 

The boundary conditions are: 1) The propellant and case remain bonded 

together, but the propellant-case assemblage is unconstrained; 2) the 
temperature of the bore surface at any instant of time is uniformly dis- 
tributed: 3) only the bore surface lo subjected to the heat flux calculated 
from Equation IV-3 (i.e., the e-.*ds of the propellant are assumed to be 
insulated); 4) the pressure at any instant of time is uniformly distributed 
over both the bore surface and the ends of the propellant. The pressure 
increases exponentially according to the relation 


P = P 

max 




(IV-10) 


Numerical Results and Discussion 


The numerical results presented here show how the propellant burning 
rate is affected by thermoelastic coupling and also how these effects are 
modified when certain parameters are varied. The results do not represent 
a detailed parametric study but are presented to show the general influence 
of certain parameters. 

For comparison purposes, a "baseline" case was established as a 
reference for making comparisons. Table IV-1 shows the basic data used 
for the baseline analysis. The results obtained for the baseline example 
are shown in Table IV-2 for 50 milliseconds. The time period chosen is 
sufficiently long to allow the dynamic effect on burning rate to decrease 
to within a few percent of the quasi-steady state rate and, as can be seen 
in Table IV-2, the influence of thermoelastic coupling on the dynamic rate 
has practically vanished. The notation for burning rate in Table IV-2 
and in later tables, r gs , rt) w / D and r t ) w represent the quasi-steady state 
rate, transient rate without thermoelastic coupling and transient rate with 
thermoelastic coupling, respectively. This table indicates that the influence 
of thermoelastic coupling is present but that its significance with regard 
to burning rate decreases very rapidly with time. 
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Table IV-1. Baseline example data. 

Burning rate coefficient (a) 

.03783 in/sec-psi n 

Burning rate exponent (n) 

.35 

Pyrolysis law (Eq. IV-4) constant (A g ) 

882.28 in/ sec 

Activation energy (E ) 

s 

28800 BTU/lb-mole 

Thermal conductivity (X) 

5.1 x 10" 6 BTU/sec- 

Density (p) 

.063657 lbm/in 3 

Specific heat a constant volume (c) 

.24 BTU/lbm-°R 

Modulus of elasticity (E) 

550 psi 

Linear coefficient of thermal expansion (a) 

5,27 x 10- 5 /°R 

Poisson's ratio (v) 

.499 

Initial pressurization rate (P^) 

14000 psi/sec 

Initial ambient pressure (P^) 

14.7 psi 

Maximum chamber pressure (P ) 
r max 

1CC0 psi 


The first parameter to be varied is the pressurization rate. This is 
done by charging the initial pressurization rate used in Eq. IV-10 from 
14,000 to 28,000 psi/sec. T*”* results obtained for this example are shown 
in Table IV-3. 

The effect on the significance of the ..nermoelastic coupling produced 
by changing the pressurization rate can be seen in Fig, IV-3. This figure 
shows the logarithm (base 10) of the absolute value of the percent difference 
between the transient rate without thermoelastic coupling and the transient 
rate with thermoelastic coupling for both the increased pressurization rate 
and baseline examples as a function of time. As can be seen, increasing 
the pressurization rate has a significant effect on the thermoelastic phenom- 
ena. This result is consistent since the thermoelastic term in Eq. IV-1 has 
a strain rate factor which is proportional to the pressurization rate. 

The second parameter altered was the elastic modulus of the propellant. 
Not only does the elastic modulus affect the strain and hence the strain 
rate, but it is also a multiplying factor in the thermoelastic term of 
Eq. IV-1. The elastic modulus was changed from 550 to 2,000 psi. This 
change is somewhat arbitrary, but it is known that a propellant does exhibit 
a considerably higher effective elastic modulus under rapid loading condi- 
tions. The results obtained for this example are shown in Table IV-4. Note 
that r ss and rf) w / 0 arc the same as in Table IV-2, since these values are 
not affected by tne elastic modulus. 

The effect of this change on the thcrmoelastic coupling is shown, 
in Fig. IV-4. which is again presented as the logarithm (base 10) 
of the absolute value of the percent difference between the 
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Table 

IV-2 . Results 

for baseline example 

Time 

r ss 

r t^w/o 

r t ) w 

(msec) 

(in/sec) 

(in/sec) 

(in/sec) 

1 

.1223 

.1122 

.1122 

2 

• 1A03 

.1A30 

.1398 

3 

.15A6 

.1820 

.1795 

A 

.1666 

.2153 

.21A0 

5 

.1770 

.2386 

.2379 

6 

.186A 

.2563 

.2559 

7 

.19A8 

.2703 

.2700 

8 

.202A 

.2813 

.2812 

9 

.2095 

.2900 

.2900 

10 

.2160 

.2969 

.2969 

11 

.2222 

.3023 

.3023 

12 

.2279 

.3065 

.3066 

13 

.2333 

,3098 

.3100 

1A 

• 238A 

.3125 

.3127 

15 

.2A33 

. 31A6 

.3156 

20 

• 26A3 

.321A 

. 322A 

25 

>2815 

.3259 

.3268 

30 

.2959 

. 330A 

.3311 

35 

.3082 

.3350 

.3357 

AO 

.3189 

.3399 

. 3A0A 

A5 

.3283 

. 3AA9 

• 3A53 

50 

.3367 

. 3AA9 

.3503 
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TaW-5 IV- 3. 

Results for modified pressurization rate example. 

Tine 

(msec) 

r 

88 

(in/sec) 

r t^w/o 

(in/sec) 

Vw 

(in/sec) 

1 

.1403 

.1122 

.1122 

2 

.1666 

.16*5 

.1510 

3 

.1864 

.2452 

.2270 

4 

.2024 

.3004 

.2946 

5 

.2160 

.3271 

.3249 

6 

.2279 

.3493 

.3479 

7 

.2384 

.3631 

.3629 

8 

.2479 

.3724 

.3729 

9 

.2565 

.3779 

.3789 

10 

.2643 

.3808 

.3822 

11 

.2716 

.3821 

.3837 

12 

.2783 

.3823 

.3840 

13 

.2846 

.3819 

.3837 

14 

.2904 

.3813 

.3830 

15 

.2959 

.3807 

.3823 

20 

.3189 

.3787 

.3800 

25 

.3367 

.3790 

. 3800 

30 

.3508 

.3810 

.3818 

35 

.3622 

.3840 

.3847 

40 

.3716 

.3875 

.3881 

45 

.3794 

.3913 

.3918 

50 

.3859 

.3952 

.3957 
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Table IV-4. Results for modified elastic modulus example 


Time 

(msec) 

r ss 

(In/ sec) 

r t>«/0 

(in/sec) 

Vw 

(In/sec) 

1 

.1223 

.1122 

.1122 

2 

.1403 

.1430 

.0930 

3 

.1546 

.1820 

.1308 

4 

.1666 

.2153 

.1777 

5 

.1770 

.2386 

.2211 

6 

.1864 

.2563 

.2424 

7 

.1948 

.2703 

.2650 

8 

.2024 

.2813 

.2785 

9 

.2095 

.2900 

.2890 

10 

.2160 

.2969 

.2972 

11 

.2222 

.3023 

.3036 

12 

.2279 

.3065 

.3086 

13 

.2333 

.3098 

.3125 

14 

.2384 

.3125 

.3156 

15 

.2433 

.3146 

.3180 

20 

.2643 

.3214 

.3249 

25 

.2815 

.3259 

.3290 

30 

.2959 

.3304 

.3329 

35 

.3082 

.3350 

.3371 

40 

.3189 

.3399 

.3416 

45 

.3283 

.3449 

.3463 

50 

.3367 

.3499 

.3510 




t w/o 


-47- 


10.00 


^ 5.00 

c 

« 

u 

M 

a 


o 

o 


x 

_ 1.00 

rH 


I 


0.50 



>4 


0.10 


0.05 


\ 



0.01 • * — • ■ 1 1 — 

0 10 20 30 40 50 

(Tine (nsec) 

Figure IV-4. Comparison of results with increased elastic modulus 
to baseline solution. 
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transient rates with and without coupling for both the baseline and 
modified elastic modulus examples. The effect of a larger modulus is 
very similar to that of an increase in pressurization rate. For this 
example , the effect J ; significantly larger than for the increased pres- 
surization rate. This should not be taken to mean that the influence of 
the elastic modulus is more significant than that of the pressurization 
rate. A much more detailed study would be required to draw any conclusion 
with regard to the relative importance of these two parameters. 

Indeed, many other parameters which affect both the burning rate and 
the thermoelastic coupling should be subjects of future parametric studies. 
Among these would be the constants in the burning rate law and the appli- 
cation of the burning rate law over a wide range of tesperatures . The 
coefficient of thermal expansion is also subject to the local temperature 
and state of strain. The experimental data required regarding the varia- 
tion in the coefficient of thermal expansion under the conditions which 
exist during the firing of a solid rocket motor has not been found. Two 
other parameters which deserve consideration are the thermal conductivity 
and Poisson's ratio. The conductivity influences the depth of the heat- 
affected zone and hence the heat transfer feedback mechanism. Poisson's 
ratio could significantly affect the state of strain if it were much less 
than O.S. 

Many of the parameters discussed above would require a considerable 
experimental effort to determine the behavior of these variables under the 
conditions experienced during the firing of a rocket motor. However, 
considerable data is available for the elastic modulus with regard to thermal 
and mechanical loading conditions. Hence, future efforts should first be 
directed toward incorporating a more rigorous model of the modulus including 
the influence of both the thermal state and loading rate. 

In sunraary, a model has been constructed to analyze the influence of 
thermoelastic coupling on propellant burning rate during transient operating 
conditions. The analysis utilizes a model of the combustion process which 
is derived from steady stare burning rate data and which is more readily 
applied than some other models considered such as presented in Ref. 19 which 
require a more detailed model of the flame structure. Examples have been 
presented which show that the thermoelnstic coupling does produce a small, 
but significant, effect during the initial portion of an ignition process. 

The results show that increasing the rate of pressurization and/or increas- 
ing the elastic modulus increases the thermoelastic effect. Several 
variables are mentioned as possible candidates for a rigorous parametric 
study and/or experimental research; in particular, the propellant modulus. 



V. PERFORMANCE PREDICTABILITY INVESTIGATION 


In designing an SRM it is important to establish the limits within 
which the performance of the SRM may be expected to be at all points in 
its operating time. These limits will be more narrow when the quality of 
the prediction is higher. Poor predictability yields broad limits which 
in turn lead to conservatism and inefficient, uneconomical designs. The 
problem addressed here is how to assess the degree of predictability so 
that limits can be set with confidence. 

The quality of a ballistic prediction is degraded by two general 
failings: 1) model inadequacies and 2) uncertainties of basic materials 

and dimensional characteristics used to evaluate the model. Examples of 
model inadequacy include inability to model exactly ignition transients 
and even equilibrium burning rate as a function of all pertinent variables. 
Likewise, for non-homogeneity of chemical species and, to a lesser extent, 
grain temperature. Even the dimensional specifications are often not 
portrayed precisely; e.g., the ovality of a nominally c.p. grain and the 
eccentricity of the grain bore with respect to the motor centerline. Al- 
though the modeling is constantly being improved, it must be admitted that 
all ballistic models leave much to be desired. 

For examples of the second general shortcoming of ballistic analysis 
which degrades predictability, uncertainties of basic material and dimen- 
sional characteristics, we have only to consider such parameters as burning 
rate, throat erosion, propellant density and again, mandrel ovality and 
misalignment. These and many other parameters are subject to processing 
and manufacturing variations. Now, although the statistical variation of 
these variables about some mean value may generally be deduced from pre- 
vious motors, for any particular design the mean value itself is a matter 
of some uncertainty until a number of motors have been built and tested. 

The uncertainty in general decreases with each additional firing. This 
explains why the predictability limits become narrower with additional 
motors tested. The problem here is to identify quantitatively the degree 
of predictability for all points in time in the development program. Per- 
haps the predictability of the first motor should be stressed; however, 
because it influences the basic design which may be impractical to change 
later. 

Application of arbitrary safety margins to allow for deviation from 
predicted performance could, of course, lead to either overly conservative 
or unsafe designs either of which is undesirable. The Monte Carlo program, 
the other hand, provides an approach for establishing limits on performance 
on a logical probability basis. 

Reference 2 shows how the Monte Carlo program can be adapted to permit 
statistical variation in mean values from motor to motor as well as statis- 
tical variation about the mean value. To evaluate the effect of the uncer- 
tainties on the limits of performance, it is only necessary to ascribe to 


on 
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each variable a variability (e.g., a standard deviation for a normal dis- 
tributed population) of its mean value corresponding to that deduced from 
experience on other programs. This variability will identify the 
uncertainty of the corresponding variable. Each variable may also be 
assigned distributions corresponding to more or less well established 
variations about the shifting (uncertain) mean values. After evaluation 
of a large number of motors selected on a probability basis in the manner 
of the Monte Carlo program, the performance limits versus time are estab- 
lished. The uncertainties can be adjusted to obtain predictability limits 
at various times in the development program as more information of the mean 
values is obtained. Also, both pair and Individual motor population per- 
formance variations may be determined. 

It is customary in ballistic predictions to adjust the constants in 
the burning rate relation used to obtain improved predictability as data 
from firings is obtained. The improvement is often dramatic. The improve- 
ment may actually be due to adjustment for inadequacies of the burning 
rate coefficient. This suggests that model inadequacies may at least 
partially be treated as uncertainties in the approach outlined here. Care 
must be taken, however, in deducing values for uncertainty from previous 
programs to segregate those that are due to actual variation in input 
parameters from those that are due to model inadequacies. 

Unfortunately, it now appears impractical to obtain the data necessary 
to establish uncertainties of any type as our queries of industry and 
government sources and survey of the literature reveals. It is most 
difficult to recover documentation on what values of a specific variable 
were used in an original design calculation versus those that were later 
established to be more valid in various stages of the development and 
production programs. This is the type of information required for the 
suggested predictability assessment. 

So our efforts to find a better way to establish predictability limits 
have proved somewhat abortive. Nevertheless, we have recorded the basic 
approach in the hope that it gives impetus to acquisition and documentation 
of more detailed statistical data on design and manufacturing variables. 



VI. CHANGES TO PREVIOUS COMPUTER PROGRAMS 


In this section, changes to the two computer programs listed in the 
Appendices to Ref. 2 are given. Each change is listed separately and 
identified with a page and line number cotinted from the top of the page. 

Brief explanations of the reason for each set of changes are also given. 

The changes are for the most part refinements in analysis or program logic. 

A few minor errors in the design analysis program have been found and are 
corrected herein. The most extensive change is the improvement of the use of 
tabular values of burning surface area during tailoff in the Monte Carlo 
program. This was accomplished earlier for the design analysis program. 

, The specific modifications are given in the form of new or revised 
cards. The cards should be punched beginning with column 7 except that state- 
ment numbers are punched to end in column 5 or unless otherwise noted. The 
card numbers given are for reference to the MSFC and Auburn University program 
listings only. Numbers for cards that have been deleted are marked with an 
asterisk. Modifications Nos. 1 through 15 apply to the Monte Carlo program 
and Nos. 16 through 22 to the design analysis program. 

Modification No. 1 - Card No. 00087050 

Purpose: To eliminate statistical analysis for single motors. 

Add card between lines 6 and 7 p. 89: 

IF(NRUNS.EQ.l) STOP 

Modification No. 2 - Card Nos. 00169800 & 00170000 

Purpose: To permit use of fractional number of slots. 

Change line 24 p. 106 to: 

501 FORMAT (6X,F10. 3, 3X,F10 .2) 

Change line 26, p. 106 beginning in card column 6 to: 

1.3,/,13X,'XTZO= \F7.3,/,13X,'S= » ,F6 .2./.13X, ’THETAG= » ,F8.5,/,13 

Modification No. 3 - Card Nos. 00077400* & 00077500* 

Purpose: To eliminate instability in mass balance iteration during 

tailoff near transition pressure (PTRAN) . 

Delete lines 6 and 7 p. 87 
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Mo dlficatlon No. 4 - Card Nos. 00069900 & 00070000 

Purpose: To eliminate occasional Instability in mass balance iteration 

caused by too small a distance burned increment at web time 

Change line 27 p. 85 to: 

IF(Y.GE.ANS.AND KOUNT.EQ.O) DELY-ANS-YLED40 .05*DELTAY 

Change line 28 p.85 to: 

IF (Y . GE . ANS . AND . KOUNT . EQ . 0 ) Y-ANS+0 . 05*DELTAY 

Modification No. 5 - Card Nos. 00004400, 00174400, 00005100, 00045800*. 

00045900*, 00046000*, 00084750 

Purpose: To improve accuracy of calculation of action time; l.e. 

time at which specified low tnrust level is reached during 
tailoff . 

Add to line 28 p. 107 and line 44 p. 71: 

,FPLOT 

Add to line 3 p. 72 
,FPLOT(999) 

Delete lines 26, 27 and 28 p.80. 

Insert between lines 31 and 32 p. 88: 

CALL INTRP1 (TPLOT, FPLOT , IPT , ATF , ATFAT ,1) 

Modification No. 6 - Card Nos. 00061950, 00061951, 00061952, 00062000 

Purpose: To eliminate unnecessary calculations when there is no 

erosive burning. 

Add 3 cards between lines 43 rad -+4 p. 83: 

3 IF (ALPHA) 131, 132, 131 
132 RN0Z=Q*PN0Z**N 
GO TO 4 

Change statement number on line 44 p. 83 from 3 to 131: 
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Modification No. 7 - Card Nos. 00002600, 00106800, 00044651, 00075000, 

00075300, 00075301, 00075900, 00075901, 00115850, 
00115851, 00115900, 00114150, 00114151, 00114200, 
00163150, 00164000, 00075600 

Purpose: To improve logic of the use of tabular values of burning 

surface area during tailoff. 

Add to lines 26 p. 71 and 18 p. 93: 

,ABTT 

Add between lines 14 and 15 p. 80: 

ABTT*=0 . 0 

Change line 30 p. 86 to: 

SUMAB=ABMA IN+ABTT 

Change line 33 p. 86 (continuation card required) to: 

SUMAB= (1 . +SUMDY / ZW-DELYW/ (2 . *ZW) ) *ABT0- (SUMDY/ZW-DELYW/ (2 .* 

ZW) ) *ABMAIN+ABTT-ABDIF1 

Change line 39 p.86 (continuation card required) to: 

SUMAB= (1 . -SUMDY/ZW+DELYW/ (2 . *ZW) ) *ABMAIN+ (SUMDY/ZW-DELYW/ (2 . *ZW) ) 
*ABTOf ABTT-ABDI FI 


Insert 2 lines between lines 12 & 13 p. 95: 


YB=Y 


IF (K.EQ.l)Y=YB-SUMDY/2. 

Prefix line 13 p. 95 with statement number: 91 

Insert 2 lines between lines 43 and 44 p. 94: 

IF0K.EQ.2) GO TO 91 
IF(K.EQ.l) Y=YB 
Change line 44 p. 94 to: 

IF(YT.I.E.Y) GO TO 8 
Insert between lines 5 and 6 p. 105: 


ABTT=ABPT+AB ST+ABNT 
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Add to lines 14 and 18 p. 105: 

-ABTT 

Change line 36 p. 86 to: 

SUMAB=ABT04-ABTT 

Modification No. 8 - Card Nos. 00081450, 00081550, 00104950, & 00104951. 

Purpose: To provide printout to indicate DP0UT and X0UT limits have 

been exceeded if appropriate. 

Insert between lines 46 and 47 p. 87: 

IF(ABS(DPCDY) .GE.DPOUT) WRITE (6,1890) 

Insert between lines 47 and 48 p. 87: 

IF (Y.GE.X0UT) WRITE(6,1891) 

Insert 2 lines between lines 41 and 42 p. 92: 

1890 FORMAT (’ DPOUT LIMIT EXCEEDED') 

1891 FORMAT (' X0UT LIMIT EXCEEDED') 

Modification No. 9 - Card Nos. 00185151, 00185152, 00191651, 00191652, 

00196550, 00196551 

Purpose: To output imbalance data to tape (see Section II). 

Insert 2 lines between lines 39 and 40 p. 109 and 8 and 9 p. Ill: 
WRITE(8, 99999) NX 

WRITE (8, 99998) (TPLOT(I) ,FDLFF(I) , IDIFF(I) , IADIFF(I) , 1=1, NX) 
Insert 2 lines between lines 9 and 10 p. 112 

99998 FORMAT (4E16 .9) 

99999 FORMAT (110) 
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Modification No. 10 - Card Nos. 00007402, 00014002, 00015202, 00015951, 

00030700, 00030800 

Purpose: To allow card input of radial grain temperature gradient. 

Insert between lines 26 and 27 p. 72: 

IREAD=3 

Insert between lines 44 and 45 p. 73 beginning in card column 1: 

C * 2 FOR CARD INPUT OF TEMPERATURE GRADIENTS * 

Insert between lines 8 and 9 p. 74: 

IF(ITEMP.EQ. 2) IREAD = 5 
Insert between lines 15 and 16 p. 74: 

IF(ITEMP.EQ. 2) ITEMP = 0 
Charge lines 19 and 20 p. 77 to: 

READ (IREAD, 3700) TUBULKO, TBULKE 

2700 READ (IREAD, 3700) (YTAB(ITAB) , TiABA(ITAB) ,TTABB(ITAB) , 
TTABC(ITAB), 

Modification No. 11 - Card Nos. 00031300, 00031500, 00032800, 00036950, 

00039500, 00077100, 00098900, 00099700, 00099701 

Purpose: To permit user to specify the extinguishment chamber pressure 

(PEXT) and the grain reference temperature (TREF) . See 
also Modification No. 13 for other necessary changes for TREF. 

Add to lines 25, 27 and 40 p. 77 and line 11 p. 79: 

, PEXT, TREF 

Insert between lines 33 and 34 p. 78: 

C * PEXT IS THE PRESSURE AT WHICH THE PROPELLANT 
EXTINGUISHES * 

Change line 3 p. 87 to: 

IF(PONOZ.LE.PEXT) GO TO 25 

Add to line 29 p. 91 inside closing parenthesis: 

,6X,F10. 2 ,6X,F10 .2 



Add to line 36 p. 91 Inside closing parenthesis (additional continuation 
card required): 

,/,13X,*PEXT- '.Fe.Z./.lSX.’TREF- * ,F6 .2 
Modification No. 12 - Card Nos. 00077000, 00077010, 00077020, 00079100 
Purpose: To improve logic of mass balance iteration during tailoff 

Change line 2 p. 87 to: 

22 IF(PBAR.LT. PTRAN) DPCDY - PBAR*DADY / (1 . -Nl) /ABAVE 
Insert 2 lines between lines 2 and 3 p. 87: 

IF (PBAR.GE. PTRAN) DPCDY - P3AR*DADY/ (1 . -N2) /ABAVE 
PONOZ *= PONO J+DPCDY *DELY 
- ^ee line 23 p. 87 to: 

A- CONTINUE 

Modification No. 13 - Card Nos. 00036960, 000A8300, 00046100, 00055000, 

00058200, 00060500, 00077600, 00048500, 00048700, 

00048800, 00048900, 00049000, 00049600, 00049700, 

00055200, 00055300, 00055400, 00055500, 00055900, 

00056000, 00058400, 00058500, 00058600, 00058700, 

00059100, 00059200, 00060700, 00060800, 00060900, 

000610CC , 00061400, 00061500, 00077800, 00077900, 

00078000, 00078100, 00078500, 00078600 

Purpose: To allow user to specify a grain reference temperature (TREF) 

See also Modification No. 11. 

Insert between lines 33 and 34 p. 78: 

C * TREF IS THE DESIGN TEMPERATURE OF THE GRAIN * 

Change oO.O in the following places to TREF: 

Line 3, p. 81; line 29 p. 80: line 22 p. 82; line 6 p. 83; line 29 
p. 83; line 8 p. 87; line 5 p. 81; lines 7, 8, 9, 10, 

16 and 17 p. 81; lines 24, 25, 26, 27, 31 and 32 p. 82; 
lines 8, 5, 10, 11, 15 and 16 p. 83; lines 31, 32, 33, 34, 

38, and 39 p. 83; lines 10, 11, 12, 13, 17 and 18 p. 87. 



-57- 


Modification No. 14 - Card Nos. 00179500, 00179520, 00179540, 00196553 

Purpose: To create a tape containing the thrust versus time data for 

each SRM. The data is written on unit 9 and the appropriate 
job control cards must be provided for this unit. This 
modification was accomplished to obtain nominal trace 
characteristics for statistical correlation of thrust im- 
balance with nominal trace characteristics (See Section II) . 

Change line 31 p. 108 to: 

2 WRITE(9, 99999) NP 

Insert 2 lines between lines 31 and 32 p. 108: 

WRITE (9, 99997) (TPLOT(I) ,FPL0T(I) ,I=1,NP) 

NP=NP+2 

Insert between lines 9 and 10 p. 112 after Modification 7 insertion: 
99997 FORMAT (2E16. 9) 

Modification No. 15 - Card Nos. 00012200, 00012500, 00012600, 00012700, 

00013900 

Purpose: To improve clarity of option descriptions. 

Insert line 26 p . 73 between OVALITY and ANALYSIS: 

OR NON-AXISYMMEl’RIC THERMAL 
Change description line 29 p. 73 to: 

FOR PLOTS, TABULAR OUTPUT AND STATISTICAL ANALYSIS 
Change description line 30 p. 73 to: 

FOR TABULAR OUTPUT AND STATISTICAL ANALYSIS 
Change description line 31 p. 73 to: 

FOR PLOTS AND STATISTICAL ANALYSIS 
Change description line 43 p. 73 to: 

FOR TAPE INPUT OF TEMPERATURE GRADIENTS 
See also Modification 10 (between lines 44 and 45 p. 73) 
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Modi fixation No. 16 - Card Nos. 00030200*. 00030300*. 00020451, 00020452, 

00020651, 00020652, 00020653, 0C020700 

Purpose: To correct error in calculation of burning rate coefficients 

Q1 and Q2 which occurs with input of grain temperature 
gradients . 

Remove Q1 and Q2 equations, lines 15 and 16 p. 158, and insert same 
between lines 13 and 14 p. 156. 

Insert 3 lines between lines 15 and 16 p. 156: 

IF(ITEMP) 10005, 10006, 10005 

10005 Q1=A1*EXP (PIPK* (1 . -Nl)* (TGR-TREF) ) 

Q?=A2*EXP (PIPK* (1 . -N2 ) *TGR-TREF) ) 

Change line 16 p. 156 to: 

10006 IF(XT.LE.0.0)GO TO 40 

Modification No. 17 - C^rd Nos. 00032151. 00032152. 00032153. 00033154. 

00Q32200 

Purpose: To eliminate iteration for burning rate wnen there is no 

erosive burning. 

Insert 4 lines between lines 34 and 35 p. 158: 

3 IF(ALPHA) 131,132,131 
132 IF(PNOZ.LE.PTRAN) RNOZ = Q1*PN0Z**N1 
IF(PNOZ.GT.PTRAN) RNOZ = Q2 ?N0Z**N2 
GO TO 4 

Change statement number on line 35 p. 158 from 3 to 131 

Modification No. 18 - Card Nos. 00037300 & 00037400 

Purpose: To eliminate occasional instability in mass balance iteration 

at web time which occurs when y (DELY) increment just prior 
to beginning of tailoff is toe small. 

Add to lines 38 and 39 p. 159: 


+0.Q5*DELTAY 



Modification So. 19 - Card No. 00045100* 


Purpose: To eliminate instability in mass balance iteration during 

tailoff when the chamber pressure crosses the transition 
pressure PTRAN. 

Delete line 20 p. 161. 

Modification No. 20 - Card Nos. 00045000, 00045010, 00045020 

Purpose: To correct error in selection of N2 and N1 when PTRAN is 

crossed during tailoff 

Change line 19 p. 161 to: 

IFCP0N0Z.LE. PTRAN) Q2=Q1 

Insert 2 lines between lines 19 and 20 p. 161: 

IF(P0N0Z.LE. PTRAN) N2=N1 
RN0Z=Q2*P0N0Z**N2 

Modification No. 21 - Card Nos. 00124800, 00124900, 00125000 

Purpose: To allow the plotting of additional data points so that the 

user may use smaller y increments. 

Change 200 to 999 throughout lines 6, 7, and 8 p. 178. 

Modification No. 22 - Card Nos. 00019810, and 0C020600 

Purpose: To correct an error in the place of calculation of the 

characteristic velocity. This error would occur when 
TRC#TREF and CSTARTyO.O. causing a small accumulative error 
in CSTAR over the operating time. 

Insert between lines 14 and 15 p. 156: 

CSTARN=CSTAR 

Change line 15 p. 156 to: 

10003 CSTAR=C$TARN*(1 .+CSTART*(TGR-TREF) ) 



VII. CONCLUDING REMARKS 


The continued research in internal ballistics performance variation 
has produced considerable improvement in the flexibility of the Monte 
Carlo and design analysis computer programs. In addition to a large number 
of minor improvements in both programs, two new computer program codes now 
permit more detailed examination of the thrust imbalance characteristics 
defined by the Monte Carlo program. Sample evaluations for the Titan IIIC 
and comparisons with flight test data show good agreement throughout the 
operating time of the motors and thus provide additional confirmation of 
the Monte Carlo method. Statistical correlation studies of the sample 
results suggest several new and possibly important approaches to predicting 
thrust imbalance. 

Design analysis program modifications now permit evaluation of the 
internal ballistic effect of general deformation of circular-perforated 
grains. Although additional validation of the model is needed, application 
of the method to two different SRMs indicates that the deformation may 
account for a substantial portion of the so-called "scale factor" on 
burning rates. 

The research in the thermoelastic coupling of propellant strain rate 
with grain temperature and burning rate has culminated in an apparently 
successful method for theoretical evaluation of the phenomenon. Assessment 
of the influence of thermoelastic coupling indicates a small but definite 
effect which warrants further investigation. 
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APPENDICES 


COMPUTER PROGRAMS 


A. Analysis of SRM Pair Imbalance Data during Specified 
Time Intervals 

B. Analysis of SKM Pair Imbalance Data at Specified Times 
during Operation 

C. SRM Design Analysis with Option Permitting Evaluation 
of Grain Deformation Effects 
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APPENDIX A 


A COMPUTER PROGRAM FOR ANALYSIS OF SRM PAIR IMBALANCE 
DATA DURING SPECIFIED TIME INTERVALS 


The computer program listed in this appendix is used to analyze a data 
tape containing motor pair data which is generated by the Monte Carlo per- 
formance analysis program described in Refs. 2 and 3. This program analyzes 
the data for a prescribed number of time intervals. The variables investi- 
gated are the absolute values of thrust imbalance and its first time deriva- 
tive, the impulse imbalance and the absolute impulse imbalance. For each 
variable the program calculates the maximum value for each pair during the 
time interval, the maximum value occurring in all time intervals for all 
pairs, the mean value for each time interval and the standard deviation for 
each time interval. The program also determines the mean value and standard 
deviation of the time at which the maximum value of the variable occurs 
within the time interval. 


Input Data 

The discussion below gives the general purpose, order and FORTRAN coding 
Information for the input data. 

Card 1 Number of time intervals to be analyzed (110) 

Col. 1-10 Number of time intervals (NSETS) 

Card 2 Time interval data (110, 2F10.4)(one card for each time interval) 
Col. 1-10 Number of motor pairs to be analyzed (NAMAX) 

11-20 Beginning time of interval (ATM1N) 

21-30 Ending time of interval (ATMAX) 

Program Listing 

The program listing is presented in Table A-l. 

Output 

The output of the program is discussed in Section II of the main report 
where a sample output illustration is also given. 
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Table A-l. Co neuter program for analysis of SRH pair imbalance data 
during specified time intervals. 

DIMENSION FD IFF ( 1000 )*TPLUT( 1000 1 ,OFOOT ( 1000 1 
DIMENSION NA.MAX160) * ATMINI SO ) * ATMAX ( 50 ) 

CALL GSIZ6I500., 11. ,11211 
READ(5,10l) NSETS 
DO 7000 IPLT=l,3 
REWIND 8 

00 6000 L*l, NSETS 

IF ( IPLT.EQ.2.AND.L.E0.1) CALL PLOT I 24 . ,0, ,-3 ) 

IF ( IPLT.E0.2.AND.L.GT.I I CALL PLOT 19. » 0. ,-3) 

IF( IPLT.EQ.3I CALL PLOT! 4. ,0. , -3) 

XMAX=0.0 
TNAX-0.0 
XO MAX =0.0 
TDMAX=0. 0 
REWIND 8 

IFIL.GT.l.ANO. IPLT.EU.l) CALL PLOT I 24. ,0 .,-3 ) 

IF( IPLT.EO.l ) READ (5.101) NAMAX ( L J , ATM IN ( L ) , ATMAX I L ) 

NM AX=NAMAX ( L ) 

T M I N = A 1 M l N ( L ) 

TMAX=ATMAX(L ) 

WR I TE ( 6 , 1 03 ) L , NMAX, TMIN,TMAX 

DO 1000 N=l,NMAX 

OF .MAX =0.0 

T0F,MAX = 0.0 

REAUI8,l0l) NX 

K0UNT=0 

K = 0 

IC=0 

DO 4000 J= 1 » NX 
CALL INPUT ( T, F , IPL T I 
IF (I .LT.TMIN) K = J 
IF(T.LT.TMIn) GO TO 4000 

1 = J-K 

IF ( T.GT. T.MAX.ANO.KUUWT . EO.O ) KUUNT=l-l 
IF ( KUUNT ) 11,11,4000 
11 TPLOT(I)=T 
FD IFF ( I )=F 
AF =A8S ( F ) 

I F ( AF.GT .UF.MAX I TDfMAX=T 

I F ( AF . GT .OF 4AX ) UFMmX = AF 

IF (DFMAX.GT.XM AX) TiMXX=TDf MAX 
IF { D r *MAX .GT . XM AX ) XMAX=DFMAX 
IF (T.LT.TMAX.ANO. J .GE.NX ) I C = l 
4000 CONTINUE 

GO TO (87,88,89), I PL T 

87 WRITE(6,I07) N , DFM AX 
GO TO 90 

88 WRITE (6*1071) N»OFMAX 


- 65 - 



Table A-l (Cont'd) 


GO TO 90 

89 WR l Tc ( b« 1072 ) N.DFMAX 

90 WR ITE(6, L06 ) TQf MAX 

CALL SIG8AR(0FMAX,Sl,S2, XMAX S, XMAXM , N , NMAX ) 

CALL SIG3AAI TDFMAX , S3 ,S4 t TMXXS , TMXXM, N , NMAX) 

IF(KOUNT.NE.O) NX=KOUNT 
IF ( IC.NE.O) NX=IC 
IF « N- 1 1 14,14,16 
14 CALL SCAlE(FDIFF,8.0»NX»1) 

IF(FDlFF(NX*i J .LT.0.0.AND.«.0*FDIFF(NX+2)+FUlFF(NX»l) .LE.O.O) 
2FlR$TV=F0lFFlNX+l )*2. 

IF(F0lFF(NX + l).LT.Q.U.AND.8.0*FDlFFiriX*2)’-FDlFF(NX-«-l).GT.0.0) 
2FIRSTV=-ABS(d.*F0IFF(NX + 2n 

IF (F01FF(NX*1 J .GE.O.O) F IRST V=-(FDIFF( NX fl I «-8.0*FDIFF (NX*2) )*2. 
DELTAV=-rIRSTV/4.0 
16 CONTINUE 

FO IFFCNX+l )=FIRSTV 
FO IFFINX+2) =JELTA\/ 

TPLCUNX*-l)=TMIN 

TPLGT(NX+ZI=(TMAX-lMlN)/5.0 

IF (N.EQ. l.ANU.L.EO.l.ANu.IPLT .EQ.Ll CaLL PLOT (4.0, 1 .5,-3) 

IF (N.NE. 1 ) GO TO 999 
GO TO C 1,2,3) , IPLT 

1 CALL AXI S( 0.0, 0.0, ’THRUST IMBALANCE ( LBF ) • ,22 , 8.0, 90.0, 
2FIASTV»0cLT AVI 

GO TO 4 

2 CALL AXISIO. 0,0.0, • IMPULSfc IMBALANCE ( LBF-ScCS ) • , 2d , 8. 0 ,90.0 , 

2F l RSTV,UELTAV) 

GO TO 4 

3 CALL AXIStJ. 0,0.0, ’ABSOLUTE IMPULSE IMBALANCE (LBF-SLCS)’, 

237,8. 0,90.0, 

2F I RS T V, UEL T AV ) 

4 CALL AXISIO. 0,4.0, ’TIME { SECS ) ' 1 1 , 5.0 , 0 .0 , TPLOT ( NX* l ) 
2 , T PLOT ( N X*2 ) ) 

999 CALL PLOT (0. ,4. ,3) 

CALL LINE< IPLUT,F01FF,Na,1,0,I) 

1000 CONTINUE 

GU TO (3,6, M , l?l I 

5 WR ITC ( 6, IOhJ XMAX 
WRITE (6, 106) TMXX 
WR ITE ( 6, 109) 

WA IT £ { 6 » 104 ) XMAXM 
WAITE(6,lll) 

WR I T E ( 6 , 104 ) XMAXS 
WR I T£ ( 6, 1 10) TMXXM 
WRITE(6,112) TMXXS 
GO TO 8 

6 WRITE16, 1041) XMAX 
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Table A-l (Cont'd) 


WR I TE ( 6 , 106 ) TMXX 
WR ITE ( 6» 109) 

WRITE (6* 1041 ) XMAXM 
WRITE<6, 111) 

WRITEIot 1041) XMAXS 
WRITE(6,ll0) TMXXM 
WR ITE (6(112) TMXXS 
GO TO 8 

7 WR 1 1 E i 6( 1 042 ) XMAX 
WR ITfc (6t 106 ) TMXX 
WRI TE (6, 109) 

WRITEI6, 1042) XMAXM 
WRITE (6, 1)1 1 
WRITE (6, 1041) XMAXS 
WRITE(6tllO) TMXXM 
WR I TE ( 6 » 1 12 ) TMXXS 
d I F ( IPLT.NE . 1 ) GU TO 6000 
REWIND 8 

DO 2000 N= l,NMAX 
UF DM AX =0.0 
T0F0MX=0.0 
RE ADI 3, 101 ) NX 
K.OUNT = 0 
K= 0 
! C =0 

DO 5000 J •= l » N X 
CALL INPUT (TfF,IPLT) 

1FIT.LT. THIN) K = J 
IF(T.LT.lrilN) GO TO 5000 
I = J-K 

IF t T . G1 .TMAX.AND.KOUNT.LU.O) KDUNT = I-1 
IF(KUUNT) 13,13,5000 
13 T PLOT I I ) =T 
FD IFF ( I i = F 

IF (T.LT.TMAX.ANU.J .GE.NX ) IC= I 
5000 CONTINUE 

1 M K DUN I .N E . 0 ) NX = KOU NT 

IF(IC.Nt.O) N X = I C 

OF DOT ( 1 ) = FO i FF ( 1 )/ 1 PLOT { 1 ) 

OF JI1AX = At) S ( DF 00 r ( 1 ) ) 

Tt)fDMX=TPLUr ( 1 ) 

DO 3000 J = 2 , NX 

OF OUT { J) = {F0IFF(J)-FDIFF( J- l ) ) / ( I PLOT ( J ) -TPLUT ( J-l) ) 
AF D= A tJ S ( OF DO T ( J ) ) 

I F ( AF 0. G T .OFDMAX ) ToFt)MX=T PL CT ( J ) 

I F ( AFD.GT .OFDMAX ) OF DMA X s A dS (OFDUT ( J ) ) 

IF (UFDMAX .GT .XOMAX ) T OMAx = T OFUMX 
3000 IF (OFDMAX. GT .XOMAX ) XUMA X=UFDMAX 
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Table A-l (Cont'd) 


MR I TE ( 6* 10B) N.OFDMAX 
W* IT£ ( 6, 1131 TUFDMX 

CALL S IGdAR ( DFOMAX , S5, S6, XDMAXS, XOMAXM ,N »NM4X) 

CALL SIGBAR( TDFDMX tS?»SB»TOMAXSt TDMAXM ,N »NMAX I 
IF IN-1 ) 15,15,17 
15 CALL SCALE ( OF DOT, 3 *0, NX , 1 ) 

IF (OF OClT (NX+li .LT . 0.0. AND. 8. 0* OF DOT (NX+2 )*DF lOTlNX+l 1 . LE.0.0 I 
2FlRSTV*0F00TCNX+l)*2. 

IF(OFDOUNXU).LT.O.U.AND.8.0*DFOOT(NX+2I+OFDOT(NXU J.GT.Q.O) 
2FIRSTV=-ABS(8.*DFDOT(NX*2) ) 

IF I OF OOT ( NX 1 1 .GE.0.0 ) F IkST V=-{DFDOT (NX + l 1*-8.0*DFDUT (NX+2) 1*2. 
OELT AV = -F IRSTV/4.0 
17 CONTINUE 

OF OUT (NX+l)=FlRSTV 

DF OUT ( NX *2 ) =DELTAV 

TPLOT(NX+l)=TMIN 

TPLOT (NX+2 ) = ( T MAX-TMIN) /20 .0 

lF(N.EQ.l.AND.IPLT.EU.l) CALL PLOT ( 9 .0 ,0 .0 ,-3 ) 

IFlN.Ew.il CALL AXISIO. 0,0.0, 'THRUST IMBALANCE RATE (L8S/SEC) ' , 31 , 
28.0,90.0,FlRSTV,DELTAV) 

IF (N. EU. DC ALL AXI S( 0.0, 4.0, ' TIME (SECS) ' ,-l l ,20. ,0.0, TPLOT (NXU) 

2 , T PLOT (NX+2 ) ) 

CALL PLOT (0.0, 4.0, 3) 

CALL L I NE ( TPLOT , OF DOT ,NX, 1,0, 1) 

2000 CONTINUE 

XOMAX 


6000 

7000 

100 

101 

102 


MR 1TE(6, 105) 

HR ITE ( 6 , 1 1 3 ) TOMAX 
WRI1E(6, 109) 
wR I TE (6 , 105 ) XOMAXM 
WRITE (6, 111) 

WR(Tb(6, 105 ) XOMAX S 
WRITE(b,ll4) TOMAXM 
W'R I T E ( 6 » 115) TDMAXS 
CONTINUE 
CONTINUE 

CALL PLCTIO. 0,0. 0.999) 
F CJ RM AT 1 4E 16 • 9 ) 

FO RM AT ( I 10, 2F10.4) 
FORMAT ( 3( 5X,E16.9) ) 


OMGWAL P AGE IS 
OF POOR QUALITY 


103 FORMAT ( 1H1 , 9X,»TMIS IS TIME INTERVAL NUMBER' ,14, //,1 OX, 'THERE ARE 
2 ',14,' SETS CiF DATA FUR THIS TIME INT Ei< VAL • , // , 10X, ' THI S TIME INT 
NERVAL BEGINS AT',F7.2,' SECS AND ENDS AT',F7.2,* SECS',//) 

104 FURMAT ( 10X, ' THE ABSOLUTE VALUE UF THE MAXIMUM THRUST IMBALANCE DUR 
2ING THIS TIME INTERVAL IS ',IHE11.4,' LBF* ) 

1041 FORMAT (10X, » THE ABSOLUTE VALUE OF THE MAXIMUM IMPULSE IMBALANCE DU 
2RING THIS TIME INTERVAL IS ',lPEll.4,‘ LBF-SECS') 

1042 FORMAT ( 10X, 'THE MAXIMUM ABSOLUTE IMPULSE 1MB 

2ALANCE DURING THIS TIME INTERVAL IS ',IP£11.4,' LBF-SECS'J 


-68- 



Table A-l (Coat'd) 


105 FORMAT! 10X* • THE ABSOLUTE VALUE OF THE MAXIMUM THRUST IMBALANCE RAT 
2E DURING THIS TIME INTERVAL IS ',IPE11.4,' LBF/SEC • ) 

106 FORMATllOX,' THIS IMBALANCE OCCURS AT * ,F7.2, * SECS',/) 

107 FORMAT (10X, • THE ABSOLUTE VALUE OF THE MAXIMUM THRUST IMBALANCE FOR 
2 MOTOR PAIR NUMBER ',14,* IS ',iPEll.4,» LBF ' ) 

1071 FORMAT !10X, ' THE ABSOLUTE VALUE OF THE MAXIMUM IMPULSE IMBALANCE FO 
2R MOTOR PAIR NUMBER ',14,' IS '.IPE11.4.' LBF-SECS') 

1072 FORMAT! 10X, 'THE MAXIMUM ABSOLUTE IMPULSE IMB 

2ALANCE FOR MOTOR PAIR NUMBER ',14,' IS ',1P£11.4,« LBF-SECS') 

108 FORMAT! 10X,' THE ABSOLUTE VALUE OF THE MAXIMUM THRUST IMBALANCE RAT 
2E FOR MOTOR PAIR NUMBER ',14,' IS ',1PE11.4,' LBF/SEC') 

109 FORMAT! 10X, • THE MEAN VALUE OF') 

110 FORMAT!/, 10X, 'THE MEAN VALUE OF THE TIME FOR THIS IMBALANCE IS ', 
2F 7.2, • SECS',/) 

111 FORMAT!/, 10X, • THE STANDARD DEVIATION OF') 

112 FORMATUOX.'THE STANDARD DEVIATION OF THE TIME FOR THIS IMBALANCE 
2 1 S ' ,F7«2 , • SECS',/) 

113 FORMATllOX,' THIS IMBALANCE RATE OCCURS AT ',F7 .2,' SECS*,/) 

114 FORMAT!/, 10X,' THE MEAN VALUE OF THE TIME FOR THIS IMBALANCE RATE l 
2S ' ,F7.2, ' SECS',/) 

115 FORMATllOX,' THE STANDARD DEVIATION OF THE TIME FOR THIS IMBALANCE 
2RATE IS *,F7.2,» SECS',/) 

STOP 

END 

SUBROUTINE SIG6AR(X,XI ,XI2,S IGX.BX, ICOUNT,N) 

XN=FLQAT { N ) 

IFIICOUNT.GT.I ) GO TO l 
XI 2=0.0 
XI =0.0 

1 XI2=XI2+X**2 
XI=XI+X 
6X=XI/XN 
XI S=XI **2 

ARG=! XI2/XN)-! XIS/XN**.") 

I F ( ARG ) 2,2,3 
c S I GX = J • 0 
GO TO 4 

3 SI GX=SUKT ( AxG) 

4 RETURN 
END 

SUBROUTINE INPUT! T,F, IPLT ) 

GO TCI (1,2,31, I PL T 

1 REA0!3,100) T,F, 1)1, D2 
GO TO 4 

2 READ! 3, 100) T,Dl,F,D2 
GO TO 4 

3 READ! 6, 100) T,Dl,D2,F 

4 RETURN 

100 FORMAT (4Elo*9) 

ENO 
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APPENDIX B 


A COMPUTER PROGRAM FOR ANALYSIS OF SRM PAIR IMBALANCE 
DATA AT SPECIFIED TIMES DURING OPERATION 


The computer program listed in this appendix is used to analyze a data 
tape containing motor pair data which is generated by the Monte Carlo 
performance analysis program described in Refs. 2 and 3. This program 
analyzes the data at a prescribed number of discrete time slices and cal- 
culates the statistical tolerance limits about a zero mean value for the 
thrust imbalance, thrust imbalance rate and impulse imbalance. For the 
absolute impulse imbalance the statistical tolerance limits are computed 
about the true mean value. 

Input Data 

The discussion below gives the general purpose, order and FORTRAN coding 
information for the input data. 

Card 1 (3F10.0) 

Col. 1-10 One sided K-f actor (SIGK1) 

11-20 TVo sided K-f actor (SIGK2) 

21-30 Initial time point to be plotted (TMIN) 

Card 2 (2110) 

Col. 1-10 Number of time slices to be taken (N) 

11-20 Number of motor pairs to be analyzed (NMAX) 

Card 3 Time Slice Card (F10.0)(one card per time slice) 

Col. 1-10 Time at which analysis is to be made (TSL) 

Program Listing 

The program listing is presented in Table B-l. 

Output 

The output of the program is discussed in Section II of the main report 
where a sample output illustration is also given. 
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Table B-l. Computer program for analysis of SRM pair imbalance data at 
specified times during operation. 

DIMENSION F{ 1000) ,FD< 1000) ,T( 1000) , TSL ( 500) , FO I S ( 500 1 , F £S (500 ) 
DIMENSION FI (500) ,FD1 (500) ,FIM(500) t FOIM{ 500) ,F ISM(50Q) 
DIMENSION S 1( 500) , S2 ( 500 ) • S3 ( 500) ,*4(500) 

CAU GSIZE(5JQ.»11«»1121 ) 

READ! 5,100) SI GK1 , SIGK2 , TMIN 
RE AD (5, 101) N, NMAX 
WRITE(6»105) NMAX,SIGK1,SIGK2,N 
10L FORMAT (2110) 

DO 4000 K* 1 » N 
4000 RE AD ( 5, 100) TSL(K) 

100 FORMAT (3F10.0 ) 

DO 99999 I PL T= l » 3 
REWIND B 

DO 3000 L=i,NMAX 
CALL luP UT ( T » F , NX , I P L T ) 

102 FORMAT (4E16. 9) 

FD(l)=F(l)/T(l) 

DO 1000 J = 2 » NX 

1000 F J(J)=!F( J)-F( J-l) )/(T(J)-T( J-i) ) 

DO 3000 K= l , N 

CALL INTRP1(F,T,NX,TSL(K),FI(K) ) 

CALL INTRPi(FD,T ,NX,TSL(K) ,FDI (K) ) 

104 FORMAT (10X,3E16*9) 

CALL SIGbARI F 1 ( K ) , Si ( K ) , S2 ( K ) , F I S ( K ) , F IM ( K ) , L , NMAX , 1PLT ) 

3000 CALL S I G b AR ( F 0 1 (K) » S 3 ( K ) ,S4 (K),FDIS(K) ,FOIM(K) ,L,NMaX, IPLT ) 

DO I 1=1, N 

FD IS( I )= SIGK2*FDIS( I ) 

IF( IPLT. ED. 3) GO TO 11 
F I S( I ) = S IGK2*F I 5 ( I ) 

GO TO l 

11 FI SMI I) = FIM( I)-SIGKl*FIS(I) 

IF ( F I SMI I ) .LE • 0*0) F I SM( I ) =0 • 0 
F I S( I ) = F I M ( I J+SlGKl^F JS( I ) 

1 CONTINUE 
DO 4 1=1, N 
WRITE (6,106) TSUI) 

IF ( IPLT .NE.i) GO TO 30 
WRITE 16,103) F IS( I ) * F (0 1 S ( I ) 

GO TO 4 

30 I F ( I PL 1 .NF..2 ) GO TO 40 
WRlTfc(6,113) FIS(I) 

GO TO 4 

40 WR I 7 1 - ( 6 » 1 3 2 ) F I M ( I ) 

WR i It (6, 123) F IS( I ) 

WRITE (6, 133) FlSM(l) 

4 CONTINUE 

IF ( IPLT— 1 ) 80,80,70 
80 CALL PL0T(4.0, L.5,-3) 
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Table B 1 (Cort'd) 


GO TO 90 

70 CALL PL0TI9. 0,0. 0,-3) 

90 CONTINUE 

CALL SCALE(Fl$,8.0,N,l) 

IF(FIS(N«-1).LT .0.0.AND.8.0*FIS(N*2)+FIS(N*1) .LE.0.0) 
2FIR$TV=FlMN*l)*2. 

IFIFISIN+I) .L T . 0.0. AND. 8.0*FIS(N+2)+FI$(N+l).GT.O.0) 
2FIRSTV=-AdS(3.*FIS(N*2) ) 

IF (FISIN+I) .Gfc.0.0) FIRSTV=-(FISlN+i)+8.0*FlS(N-*-2) )*2. 

0ELTAV=-FlRSTV/4.0 

FISIN+1) =FIRSTV 

F I S( N+2) =0ELTA V 

TSLIN + l ) =TMIN 

TSL(N+2)=(TSL(N)-TSL(N+1) )/5.0 
IF I IPLT.NE • L ) GO TO 10 

CALL AXI S( 0.0, 0.0, ' THRUST IMBALANCE I LBF ) • ,22 , 8.0 ,90.0, 
2FIRSTV,DELTAV) 

GO TO 15 

10 I F ( IPLT.NE. 2) GU TU 20 

CALL AXISIO.0,0.0, • IMPULSE IMBALANCE I L BF- StCS ) ' , 28, J. 0,90.0, 

2FI RSTV.OfcL TAV) 

GO TO 15 

20 CALL AXISIO. 0,0.0, 'ABSOLUTE IMPULSE IMBALANCE ILBF-SECS I » . 37, 8.0, 

290. 0, F IRSTVtOELTAV) 

15 CALL AXI SIO. 0,4.0, • TIME (SECS) * ,-11,5.0,0.0, TSLIN+l) 

2 , T SL I N+2 ) ) 

CALL PLOT (0., 4. ,3) 

CALL LINEITSL, FIS, N, 1,0,1) 

DO 2 I = l , N 
FI St I )=-FISCI) 

2 IF; IPLT.EQ.3) FISI I ) =F I SMt I) 

CALL PLOTIO. ,4. ,3) 

CALL LINEITSL, F IS,, 0,1 ) 

IF( 1PLT.NE.1) GO T. 999 

CALL PLOT 19.0,0.0, • 

CALL SCALE! F0IS,8. , ,1) 

»F (FOISI N+l ) .LT.O.v ,AND.8.0*F0IS I N+?) + FU ISIN + l ) .LE.0.0) 
2HRSTV»FOIS(N+l)*2. 

I F (FOIS (N + l) .LT .0.0.ANu.6.0*f 0 l S (N+2) + F0 I S I N ♦ 1 ) .GT.U.O) 
2FUSTV»-Ai)S(3.*F0IS(N+2) ) 

IF (FI) IS (N + l) .Gt.O.O) FIRST V= - 1 F 0 I SIN+1J+ 8. 0*FD I SI N+2) )*2. 

DELTAV=-FIRSTV/4.0 

F0IS(N+1)=FIRSTV 

FOI$(N+2)=OtLTAV 

CALL AXI SI 0.0, 0.0, ' THRUST IMBALANCE RATE I LBF/ SEC ) 1 , 31 , 

28. 0. 90. 0, F I RSTV,DELTAV) 

CALL AXISIO. 0,4.0, 'TIME I SEC S ) • , - 1 1 , 5 . ,0 .0 , T SL I N + 1 ) 

2 , T SL ( N+2 ) ) 


- 72 - 



Table B-l (Coat'd) 


CALL PL0T(0.0,4.0,3) 

CALL LINE(TSL,FDIS«N, 1,0,1) 

00 3 1=1, N 
3 FOISl 1)*-FDISI l) 

CALL PLOT (0.0, 4 .0,3) 

CALL LINEITSL.FOIS.N. 1,0,11 
99999 IFC IPLT.EQ.3I CALL PLOTCO. 0,0. 0,999) 

RETURN 

103 FORMAT ( 10X, 

2* THE ♦ OR - K2* SIGMA LIMIT ABOUT A ZERO MEAN FOR THE THRUST IMBALA 
2NCE IS * , 1PE 11.4, • LBF*»/,10X,'THE ♦ OR - K2*SIGMA LIMIT ABOUT A Z 
?ERO MEAN FOR THE THRUST IMBALANCE RATE IS '.IPE11.4,' L8F/SEC * »// ) 

105 F0RMATCIH1.9X, 'THERE ARE *,I4,« SETS OF MOTOR PAIR DATA* ,/, 10X, 1 TH 
2E ONE SIDED K FACTOR CK1 ) FOR THIS SAMPLE SIZE IS «,F 7.3,/, 

313X, • THE T «tO SIDED K FACTOR IK2) FOR THIS SAMPLE SIZE IS ',F7.3,/ 
2, 1 OX, 'THE RESULTS ARE CALCULATED AT ',14,' TIME SLICES',///) 

106 FORMAT ( I OX, • TM I S TIME SLICE HAS TAKEN AT »,F7.2,' SECS' I 
113 FORMATCIOX, 

2 'THE ♦ OR - K2*SIGMA LMIT ABOUT A ZERO MEAN FOR THE IMPULSE IMBAL 
2ANCE IS ',iPE11.4,' LBF SECS',//) 

123 FORMAT ( 10X, 

2 'THE ♦ K1*S I GMA LIMIT ABOUT THE MEAN FOR THE ABSOLUTE IMPULSE 
2 IMBALANCE IS ',1PE11.4,' LBF- SEC * ) 

133 format; 10X, 

2 'THE - <1*SIGMA LIMIT ABOUT THE MEAN FUR THE ABSOLUTE IMPULSE 

2 IMBALANCE IS '.lPEll.A,' LBF-SEC * » // ) 

132 FORMAT ( 1UX, * THE MEAN VALLE OF THE ABSOLUTE IMPULSE IMBALANCE IS ', 
21PE 1 1.4,* LBF -SECS') 

STOP 

END 

SUBROUI INE SI GBAR ( X , X I ,X 12 , S 1GX, BX, 1C0UNT ,N, IPLT) 

XN=FLOAT IN) 

IF ( I COUNT . GT * I ) GO TO 1 
XI 2=0 .0 
XI =0.0 

1 XI 2=XI2«-X**2 
Xt=XI+X 

BX =X I /XN 
Xl S=X1**2 
IFIIPLT-3) 2,3,2 

2 ARG= ( XI 2/XN ) 

so ;o 4 

3 ARG=(XI2/XN)-( XlS/XN**Z) 
tr (ARG.LE.O.QJ ARG=0.0 

4 SI GX=SJRT ( ARG) 

RETURN 

END 
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Table B-l (Cont'd) 


SUBROUTINE INTRPl C Y , T.N, 1 T ,DY I 
01 MENS IUN YIN), TIN) 

Nl=N-l 
OY =0 • 0 
00 l 1*1, Nl 

IFITT.GE.TIII.ANO.TT.LT.TI !♦!)) 0Y= U Y 1 1 ♦ !)- YII) ) /I Tl I ♦ll-Tl I) )) 
2*1 TT-T 1 1 ) ) +Y( I ) 

IFIDY.NE.O.O- RETURN 
1 CONTINUE 
return 

ENO 

SUBROUTINE INPUT IT , OX, NX, I PLT) 

DIMENSION TIIOOO) ,0X11000) 

READ! 8,100) NX 
GO TO (1,2,3) , IPLT 

1 RE ADI 8,200) I Ti I ) ,0X1 1 ) ,01 ,02 , 1 = 1 ,NX) 

GO TO <i 

2 RE AUl 8,200) ( T ( I ) ,01 ,0X1 I ) ,02 » I =1 ,NX ) 

GO TO 4 

3 RE ADI 8,200) I T 1 1 ) , 01 ,02 , OX I I ) , 1 = 1 ,NX ) 

4 RETURN 

100 FORMAT 1 1 10) 

200 FORMAT I 4E 16. 9) 

ENO 
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APPENDIX C 


THE SRM DESIGN ANALYSIS PROGRAM WITH OPTION PERMITTING 
EVALUATION OF GRAIN DEFORMATION EFFECTS 


This appendix contains the instructions for the preparation and 
arrangement of the data cards. In addition, a conf>lete listing of the 
program statements is given. The program was written for use on an IBM 
370/155 computer and requires approximately 100K storage locations on that 
machine. The program also is designed to be used with a CALCOMP 663 
drum plotter. The plotter requires one external storage device (magnetic 
tape or disk). However, only minor program modifications are required to 
eliminate the plotting capability of the program. 


Input Data 


The discussion below gives the general purpose, order and FORTRAN 
coding information for the input data. All of the Xu.ta, except for the 
changes described herein, are identical to the input data d" scribed in 
Appendix B of Ref. 2. 

Card 4A Replaces Card 4A of Ref. 2. Ignition, inert weight and grain 
deformation options (4X,I1,9X,I1,9X,I1) 


Col. 


1-4 IGO = 


f 0 

For 

5 f 


l 1 

For 

6-14 I WO 

= 

fo 

For 

15 ) 


U 

For 

16-24 ISO 

- 

f° 

For 

25 1 


U 

For 


no ignition calculations 
ignition calculations 

no inert weight calculations 
inert weight calculations 

no grain deformation effects 
grain deformation effects 


Card 9 Replaces Card 9 of Ref. 2. Uniform temperature and giain defor- 
mation constants (3 cards). 


Card 9A Uniform temperature card (input only if ITBMP=1) (5X.F10.0) 


Col. 1-5 TGR = 


6-15 Value of TGR 
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Card 9B Grain deformation constants (input only if ISO * 1) (7X,E12.4, 


10X,E12.4,9X,F8.5,9X,F8.5) 


1-7 

PMOD - 


8-19 

Value of 

PMOD 

20-29 

CMOD - 


30-41 

Value of CMOD 

42-50 

PMU « 


51-58 

Value of 

PMU 

59-67 

CMU - 


68-75 

Value of 

CMU 


Card 9C Grain deformation constants (input only if IStW.) (10X.E12.4, 


10X,F8.5) 


Col. 1-10 

ALPTS = 

11-22 

Value of ALPTS 

23-32 

TALC = 

33-40 

Value of TAUC 


Table C-l represents an example set of data. Table C-2 is a sample 
computer printout obtained with this input data. 

Program Listing 

Table C-3 presents the complete program listing. The program has been 
designed to produce graphical representations of the computatior.al results 
on the CALCOMP plotter. To delete the plotter compilation requirements, 
dummy subroutines may be substituted for the following subroutines: 

CSIZE, PLOT, LINE, AXIS, and SYMBOL. 
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Table C-l. Example data sheets for design analysis program 


with grain deformation (Cont’d). 








Table C-2. Sample computer printout for design analysis program with grain deformation. 


IABLLAR VALLtS FCR VI 6CLAL ll » C A6*C IN 
a-rPA. 0 . C ABSK* 0.0 


ABNK- 1 .27376 03 


APHK. 0.0 


APNK« 0.0 


VC I T > 1.9CCCE 03 


• ECUlLICRlLiP BLRMNG •♦•* 


INITIAL ..EYNCLCS Nli^Bte* 9.2779E 07 


T 1 vi * 

c.c 


Y- 

0.0 

S'. 01* 

2.A7ACE- 

■31 

atEAS* 

i. 7U1E-C1 

PTM* 

t .C 355c 

C .a 

i\ c z * 

6.61S3E-SI 

PAT*. 

I.AC340 

Cl 

C F V A C * 

1.72E66 QC 

TS?* 

2. IE 75E 

C2 

C F » 

I.A333E 00 

CfvO* 

l.t«59E 

0 0 

ItCT* 

C.C 

*.?• 

• A 


RAC * R * 

1.2C93E-37 

cr. 

8.250.6 

LG 

AiM t-6.* 

A.A179E Cl 

CEL* 

l . 3 EC 5 E 

CC 



IThtlA 

• A.A762E- 

02 KhC5* 2.C.66E- 


P C N C 2 • 

3.813AS 

32 

Pf-EAC* 

6.6915* 

32 

SLVAE* 

7.38A30 

03 

SC* 

0.3 


FVAC* 

J.3B76E 

CA 

F« 

2.7 T«**-e 

06 

VC* 

1.0827E 

GA 

»CCT* 

1.2706E 

02 

IIVAC* 

0.3 


ISPVAC* 

2.66586 

02 

EPS* 

7.6EA0E 

00 

ALT* 

O.C 


APNC2 • 

5.59A7E 

Cl 

CCF* 

1.977.1* 

Ov 


J RA7P* 0.18716-01 RATR2* 5.9L3SE-02 


IETH> 1.000GE 00 V*TA« 0.0 


1ASLLAR VALLES FC« VT» l.OCC *1*0 I* 

A?m. 0.0 ABSK* 0.0 A8M* 1.32736 03 **Ht» 0.0 


APNKa 0.0 


TIP:* 

AC. 17 


V* 11.63 


9-.C2* 

c.C 


RhEAG* 

0.3 


Pfi** 

1.2803E 

Cl 

P\C1- 

A.tilAE 

-C2 

PAT». 

1 . A t 5 4c 

01 

CF VAC • 

1.7227E 

CC 

ISP. 

o .c 


CM 

•1 . IAcBE 

OA 

crvo* 

t .tjtrs 

c c 

HOT* 

1.94546 

J4 

»p* 

2.2200c 

C 3 

3ALEA* 

3. 7SjIc 

-03 

c: * 

d.*647: 

:: 

iPPcAC* 

7.2E30E 

02 

C T C * 

- 1 . 1 A83E 

CA 





E T ► £ ! A* 2..5A5E-03 *hC5« 2.C17SE- 


*P 1 * B. 220 CE 03 

wpj- a.ji'.tt 03 

up * 3 . i ^ C _ E 03 

PhVAA. t . S 7 5 It 02 
I SP* l.ltili C 2 
101 'vaC* 2 . 451 A 6 0 2 

| I C T * t. 9 AA 6 c 06 
I T VAC * i.liO’E 06 
f Av» A. 8 AC 70 OA 
fVAOAV. 5 .A 374 E CA 
PL*. A *» ■ J ,02 i 4 c C 2 
VC 1 * 1 . 0 62 7 0 GA 
VCF* 1 . 3 7 FOE 05 
LAPBCA* 9 . 21 A 3 E -01 


PCNC2* 

9.1 37c E-0 3 

PP6A0* 

9.137 E-OI 

slvab* 

8.95A5E C 3 

SC* 

0.3 

F VAC • 

e.6599E-01 

F* 

0.0 

VC- 

1.3780E C5 

peer* 

3.266 JE-C1 

ITVAC* 

2.UA3E 36 

ISPVAC* 

2.651 E 02 

EPS* 

7.1A36E CC 

ALT* 

0.3 

APNC2* 

7. A119E 02 

CCF* 

1.56316 CO 

1 RATP 

• 9.99976- 

01 RATR2 

• 9.91786" 


XETHa l.OCOOfc 00 VET*» I.IJT6E 01 
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TABLE C-3 


******************************************************************* 

SRM DESIGN AND PERFORMANCE ANALYSIS 
PREPAREO AT AUBURN UNIVERSITY 
UNDER NOD. NO. 14 TO COOPERATIVE AGREEMENT MITH 
NASA MARSHALL SPACE FLIGHT CENTER 

BY 

R. H. SFORZINI ANO M. A. FOSTER* JR. 

AEROSPACE ENGINEERING DEPARTMENT 
SEPTEMBER 1975 


MODIFIED FOR GRAIN OEFORMAT ION INTERNAL BALLISTIC EFFECTS 

» 

• BY 

R. H. SFORZINI AND DAVIO F. SMITH 
i AEROSPACE ENGINEERING DEPARTMENT 

• SEPTEMBER 1976 

t**************************** *********************** **************** 


INTEGER GRAIN 

REAL MGEN*MDIS * MNOZtMNl *JROCK*N*L*MEl*ME*ISP*ITOT * MU* MASS* I SPVAC 
REAL N1 «N2, NSEG ,K1 *K2 * KEH* KEN* NS *LCC * LTAP 
REAL M2* MDBAR* I SP2. I TVAC ,KA, KB ,L AMBDA , I T V 
COMMON/CONST 1/ZW* AE » AT , THE TA , ALF AN 

COMMON/CONS T2/CAPGAM , ME * BOTE , ZE T AF * TB , HB .GAME , CGAME , TOPE , ZAPE 

COMMON/CONST3/S,NS.GRAIN,NTABY,NCARD 

COMMON/CGNST4/DELOI *0O»ZO*Dl 

COMMON/ VAR I Ai/Y* T » DELY* DELTA T * PONOZ * PHEAD»RNOZ »RHEAO * SUMAB* PHMAX 
COMMON/ VARI A2/ABPORT, A8SLOT , ABNOZ , APHE AO , APNOZ , DAUY , ABP2, ABNi , ABS2 
COMMON/VAR I A3/ I TOT , I T VAC * JROCK , I SP , I S P V AC , MO I S , MNOZ , SG , iUMKT 
COMMON/ VAR I A4/RNT »KHT,SUM2»Kl»R2,R3*RHAVE»f<NAVE » kBAR » YB * KOUNT » TL 
COMMON/VAR I A5/ABMA IN, AB TO , SUMDY , VC I » Ad T T ,PT RAN 

COMMON/VARIA6/WP2,CF,WP,RAOER,EPS, VC.FLAST .TLAST, DT, PONTOT , WP l 
COMMON/ VARI A7/TIMX,FV, I TV, NX 
COMMON/VARI A8/YDI 

* COMMON/VARI A9/EIHETA* RHG5 *RAT P , PMCD * CMCD • PMU * CMU* AL PT S • RATR2 

► COMMON /VARI20/Yt'TA»XETH* ISO 

COMMON/IGNl/KA,KH,UFS,RHO,L,PMIG,TI 1 » T 12 ,CSIG,01*NI*02*N2 
CUMMON/IGN2/ALPHA,dfcIA,PBIG,kklG,CfcLT JG,X*TOP, ZAP 
CQMMON/PLOTT /NUMFL T t 161 * IPO* NUUM, IPT.IGP 
DIMENSION YTABI30) ,TTA8(30) 

DATA PI,G/3.i4l59,32.i ? 25/ 

CALL GSIZE (416., 1 l.D, 1100) 

CALL PL0T16.25,2. ,-3) 

»OP=0 

RE ADI 5,500) NRUNS 
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Table C-3 (Cont'd) 


C *****************************'*************************************** 
C * READ IN THE NUMBER OF CONFIGURATIONS TO BE TESTED * 
C ********************************************************************* 


1902 

1901 


NTABY=0 

NCARD-0 

00 901 I=1»NRUNS 

NEXTR=NTABY-NCARD 

IF I NEXTR 1 1901, 1901 ,1902 

READ! 5, 1903) (01, 02,03,09,05, 06, 1 EX-1, NEXTR) 
WRITE (6,602) 1 
RE A0(5, 11111) NT AS , NTABY 
RE A0( 5,999) SUMOY , ANS s ZH» Y * T , DELI l 

lnT.CHT.RUT.BI.RP.U 4.RHA1/P.BMAUF.MI 


r AT . RMfO’ .JO! 


******************************************************************** 


SET INITIAL VALUES OF SELECTED VARIABLES EQUAL TO ZERO * 
♦♦♦NOTE*** THESE VALUES MUST BE ZEROED AT THE BEGINNING OF * 
EACH CONFIGURATION RUN * 

******************************************************************** 


READ (5,991) IG0,IW0,IS0 

READ! 5,993) IPO, ( NUMPLTt J J) , JJ = l , 16 1 , I TEMP 

******************************************************************** 


READ IN THE USER'S OPTIONS ♦ 

* 

VALUES FOR I GO ARE ♦ 

0 FOR NO IGNITION TRANSIENT CALCULATIONS ♦ 

1 FOR IGNITION TRANSIENT CALCULATIONS * 

VALUES FOR I WO ARE ♦ 

0 FOR NO INERT WEIGHT CALCULATIONS ♦ 

1 FOR INERT WEIGHT CALCULATIONS * 

VALUES FOR ISO ARE ♦ 

0 FOR NO GRAIN DEFORMATION EFFECTS * 

1 FOR GRAIN DEFORMATION BALLISTIC EFFECTS * 

♦♦♦NOTE*** IF GRAI N=2 ( ALL STAR GRAIN), ISO MUST BE 0 * 

VALUES FOR IPO ARE * 

0 FOR NO PLOTS * 

1 FOR PLOTS OF EQUILIBRIUM BURNING ONLY * 

2 FOR PLOTS OF IGNITIGo. TRANSIENT ONLY * 

3 FOR PLOTS OF B'lTH IGNITIUN TRANSIENT AND * 

EQUILIBRIUM BURLING * 

VALUES FUk NUMPLT(JJ> Ak_ (NOT REQUIRE!) FCP. 1P0 = U) * 

0 IF SPECIFIC PLOT IS NUT DESIRED ♦ 

1 *F SPECIFIC PLOT IS DESIRED * 

ORDER OF SPECIFICATION OF NUMPLTIJJ) IS * 

1 PHEAD VF TIME * 

2 PONOZ VS TIME * 

3 PHEAD AND rO* OZ VS TIME * 

9 RHEAD VS TIf- * 

5 RNOZ VS TIME ♦ 
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Table C-3 (Cont’d) 


C * 6 RHEAD ANO RNCZ VS TIME * 

1000 CONTINUE 

7 SUMAB VS TIME * 

6 SG VS TIME * 

9 SUMAB ANO SG VS TIME * 

10 F VS TIME * 

11 FVAC VS TIME ♦ 

12 F ANO FVAC VS TIME * 

13 VC VS TIME * 

14 SUMAB VS YB ♦ 

15 SG VS YB * 

16 SUMAB ANO SG VS YB * 

VALUES FOR ITEMP ARE * 

0 FOR TEMPERATURE GRAOIENT * 

1 FOR UNIFORM TEMPERATURE * 

NTAB IS THE NUMBER OF Y STATIONS FOR WHICH TABULAR * 

TEMPERATURES ARE SPECIFIED * 

NTABY IS THE NUMBER OF Y STATIONS FOR WHICH TABULAR AREAS * 
ARE SPECIFIED * 

************* ******* ******** **************************************** 
► WRITE 16,492) I GO, I WO, I SO 

WRITE (6,494) I PO, ( NUMPl T ( J J ) , J J=l, 161 , ITEMP 
WRITE(6, 11112) NTAB, NTABY 

RE ADI 5,501) RN2Ni,RH0,Al,Nl, ALPHA , BET A ,MU ,CSTAR 
C ********************************************************* ************ 


C * RE AO IN BASIC PROPELLANT CHARACTERISTICS * 
C * ♦ 
C**»* MODIFICATION MADE 5/3/76 

C * RN2NI IS THE RATIO OF THE NOMINAL VALUES OF THE BURNING RATE * 
C * EXPONENTS ABOVE AND BELOW THE TRANSITION PRESSURE * 
C * RHO IS THE DENSITY OF THE PROPELLANT IN LBM/IN**3 * 
C * A1 IS THE BURNING RATE COEFFICIENT BELOW THE TRANSITIUN * 
C * PRESSURE * 
C ♦ N1 IS THE BURNING RATE EXPONENT 8ELCW THE TRANSITION PRESSURE * 
C * ALPHA ANO BETA ARE THE CONSTANTS IN THE EROSIVE BURNING * 
C * RELATION OF ROBILLARC ANO LENOIR * 
C * MU IS THE VISCOSITY OF THE PROPEL! ANT GASES * 
C * CSTAR IS THE CHARACTER 1ST IC EXHAUST VELOCITY IN FT/SfC * 


C ** ********************** *************************** ****************** 
WRITE J 6, 603) RHC,Al,Nl,ALPHA,BETA,f'U,CSTAR,RN2Nl 

KH0=RH0/32. 174 

READ! 5,502) L , TAU, DE , DT I , THE T A , ALF AN, L TAP, XT , ZO, C ST ART , PTRAN 
C ********************************************************************* 


~ * READ IN BASIC MOTOR DIMENSIONS * 
C * * 
C * L IS THE TOTAL LENGTH OF THE GRAIN IN INCHES * 
C * TAU IS THE AVERAGE WEB THICKNESS OF THE CONTROLLING GRAIN * 
° * LENGTH IN INCHES * 
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Table C-3 (Cont'd) 


OE IS THE DIAMETER OF THE NOZZLE EXIT 1.4 INCHES * 

OTI IS THE INITIAL DIAMETER OF THE NOZZLE THROAT IN INCHES * 
THETA IS THE CANT ANGLE OF THE NOZZLE WITH RESPECT TO THE ♦ 
MOTOR AXIS IN DEGREES * 

ALFAN IS THE EXIT HALF ANGLE OF THE NCZZLE IN DEGREES * 

LTAP IS THE LENGTH OF THE GRAIN AT THE NOZZLE END HAVING * 

ADDITIONAL TAPER NOT REPRESENTED BY ZC IN INCHES * 

XT IS THE DIFFERENCE IN WEB THICKNESS ASSOCIATED WITH LTAP * 

ZO IS THE INITIAL DIFFERENCE BETWEEN WEB THICKNESSES AT THE ♦ 

HEAD AND AFT ENDS OF THE CONTROLLING GRAIN LENGTH * 

OSTART IS THE TEMPERATURE SENSITIVITY OF CSTAR * 

AT CONSTANT PRESSURE * 

PTRAN IS THE PRESSURE ABOVE WHICH THE BURNING RATE EXPONENT * 
CHANGES ♦ 


******************************************************************** 
N2-Ni*RN2Nl 
A2=A1*PTRAN**( NI-N2) 

WRITE (6, 604) L *TAU»D£»DTI , TH ETA , ALFAN , LT AP , XT , ZO,CST ART , PTRAN ,N2 
THETA=THETA/57.2V57B 
AL FAN=ALFAN/ 57 .29578 

RE A0( 5*503 ) DELTAY , XOUT , DPOUT ,ZETAF ,T B ,HB , GAM.ERREF , PREF , 
1DTREF,PIPK,TREF,GAME,PLXT 
IF{ ITEMP.NE.O) GO TO 10000 

RE ADI 5, 700) ( YTA6I I TAB) » TT AB I I TAB) , I T AB= 1 , NT AB ) 

WRITE I 6, 701) 4Y1ABI ITAB) ,TTAB( ITAB) »ITAB=1»NTA8) 

GO TO 10004 
10000 RE ADI 5 « 10001 1 TGR 

C ****************************************** *************************** 


C * READ IN BASIC PERFORMANCE CONSTANTS * 
C * * 
C * DELTAY IS THE DESIREO BURN INCREMENT GURING TAILGFF IN INCHES * 
C * XOUT IS THE DISTANCE BURNEO IN INCHES AT WHICH THE PRUPELLANT * 
C * BREAKS UP ♦ 
C * OPOUT IS 'HE DEPRESSURIZATION RATE IN LB/IN**3 AT WHICH THE * 
C * PROPELLANT IS EXTINGUISHED * 
C * ZETAF IS THE THRUST LOSS COEFFICIENT * 
C * T3 IS THE ESTIMATED BURN TIME IN SECONDS * 
C * HB IS THE ESTIMATED BURNOUT ALTITUDE IN FEET * 
C * A2 IS THE BURNING PATE COEFFICIENT ADCVE THE TRANSITION * 
C * PRESSURE * 
C * GAM IS THE KATIG OF SPECIFIC HEATS FOR THE PROPELLENT GASES * 
C * EKREF IS THE REFERENCE IHRUAT EROSION KATE IN IN/SEt ♦ 
C * TGR IS THE TEMPERATURE OE THE GRAIN IN DEGREES F * 
C * PREF IS THE REFERENCE NUZZLt STAGNATICN PRESSURE IN LB/IN**2 * 


C * OTREF IS THE REFERENCE THROAT DIAMETER IN INCHES * 
C * PIPK IS I HE TEMPERATURE SENSITIVITY COEFFICIENT UF PRESSURE * 
C * AT CONSTANT K PER DEGREE F * 


C * TREF IS THE OES IGN TEMPERATURE UF THE GRAIN IN DEGREES F * 


- 83 - 



Table C-3 (Cont’d) 


C * GAME IS THE EFFECTIVE GAMMA AT THE NOZZLE EXIT PLANE ♦ 

C ♦ PEXT IS THE PRESSURE AT WHICH THE PROPELLANT EXTINGUISHES IN * 

C * LB/IN**2 * 

- C ************* 4 ******************************************************* 
10004 WRITE! 6* 606 1 DELTA Y,X0UT *OPGUT* ZETAF *TB «HB * GAM* £RREF*PREF» DTREF 
1 * P IPX* A2 * TREF* GAME * PEXT 
IF | ITEMP.NE .0) WRITEI6, 100021 TGR 
NCARD-0 
N0UM=0 
IPT-0 
MN la, 85 
Z=ZQ 
S*0.0 
NS*0,0 
KO UNT=0 
ABMAIN=0.0 
ABT0=0.0 
ABTT*0, 

TL AST= l « 

DELY=OELTAY 

► ET HETA=0. 

* PMU=U* 5 

* ALPTS=0. 

► Y E T A= 0 « 

► 0YETA=0, 

► PMOD= 10000, 

► YE =0,0 

► XETH=0. 

TQP=GAM+ 1, 

BOT=GAM-I. 

ZAP=TUP/! 2,*B0T) 

CAPGAM=S(3RT 1 GAM )*( 2. /TUP >**ZAP 
V0PE=GAME+1. 

B0TE=GAME-1, 

ZAPE=TOPE/C 2,* BOTE ) 

CGAME= SORT! GAME) *12. /TOPE )**ZAPE 
AE=PI*0E*0E/4. 

CS T ARN=CS T AR 
1 IF(XT.LE.O.O) U=0.0 

r .f I ITEMP.NE.O) GO TO 10003 

CALL INTRPlITTA8,YTAd»NTAB»Y*TGR,0) 

•;***■* MODIFICATION MADE 1/7/76 

Ql=Al*EXP(PIPK*l 1.-N11*! TGR-TREF) 1 
Q2=A2*EXP! PIPK*! 1.-N2)*! TGR-TREF ) ) 

WRITEI6* /Ol) Y » TGR 

10003 CS T AR=CST ARN* ( l.*C START*! TGR-TREF I ) 

C**** MODIFICATION MADE 1/7/76 

IF! IT EMP 110005*10006, 10005 
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Table C-3 <'Cont'd) 


10005 Qi-Al*EXP(PIPK*( l.-Ni)*( TGR-TREF) ) 

Q2=A2*EXPIPIPK*l l.-N?>*(TGR-TREF) I 

10006 IF (XT. LE. 0.0) GO TO 40 
TL S ( Y— TAU+XT +Z/2* ) *L TAP /XT 
IF <TL«LE.O.O) TL=0.0 
1FITL.GE.LTAP) TL=LTAP 

40 IF ( T I 41,41,42 

41 OT=DT I 
GO TO 43 

42 RAOER=ERREF*( ( PONOZ/PREF ) **0 .8 ) *t IDTREF/DT)**0. 21 
DT=DT«‘(?.0*RADER*DELTAT) 

43 AT=PI*0T*0T/4. 

EP S=AE/ AT 

IFC IGO.EQ.O.OR.Y.GT.O.O) GO TO 900 

READ! 5, 97) K A, KB, UFS ,CS IG - PM 1G , T I I , T I 2 ,Rft I G , CELT IG, PBI G 

C ******** ************ ****>•-* ******** ****** ******* ************* ********* 


C * READ IN VALUES REQUIRED FOR IGNITION CALCULATIONS * 
C * ***NOT E*** NOT REQUIRED IF IG0=0 * 
C * * 
C * KA AND KB DEFINE THE CHARACTERISTIC VELOCITY IN FT/SEC * 
C * CSTR = KA + KB * PRESSURE * 
C * UFS IS THE FLAME-SPREADING SPEED IN IN/SEC * 
C * CSIG IS THE CHARACTERISTIC VELOCITY OF THE IGNITER IN FT/SEC * 
C * PM IG IS THE MAXIMUM IGNITER PRESSURE IN LBS/IN**2 * 
C * Til IS THE TIME OF MAXIMUM IGNITER PRESSURE IN SECONDS * 
C * T 1 2 IS THE TIME! IN SECONDS) FOR ’HE IGNITER PRESSURE TO * 
C * DROP TO 10 PER CENT OF MAXIMUM VALUE(PMIG) * 
C * RRIG IS THE AVERAGE REGRESSION RATE OF THE FIRST HALF OF THE * 
C * IGNITER PRESSURE TIME TRACE IN LBS/I N**2/SEC * 
C * DELTIG IS THE TIME INCREMENT FOR IGNITION TRANSIENT * 
C * CALCULATIONS IN SECONDS * 
C * PBIG IS THE BLOWOUT PRESSURE OF THE MAIN MOTOR BLOWOUT PLUG * 
C * IN LBS/IN**2 * 


q ********************************************************************* 
WR ITE 16,842) K A , KB , UF S ,C S I G , PM I G , T I 1 , T I 2 ,RR I G, CEL T I G , PB I G 
900 IF ( IWQ.EQ.O.OK.Y.GT.O.O) GO TO Bi2 

RE ADI 5,600) OTEMP, S I GMAP , S I GMAS , X 1 , X2 , S YCNGM , DCC , P S I C , DELC , LC 

lC,NSEG,HCN,SYNNOM, PS IS, PS I A,K1 ,K2,PSI I NS , 0 1 L i NS , K EH , KEN , US INEK , T AU 
2L, WA 

C **********************************^* ********* ************************ 


C * READ IN BASIC PROPERTIES REQUIRED FOR WEIGHT CALCULATIONS * 

* ***NOTE*** NOT REQUIRED IF IWO=0 * 

c * * 

c * DTEMP IS THF MAX EXPECTED INCREASE IN TEMPERATURE ABOVE * 

C * CONOITIUNS UNDER WHICH MAIN TRACE WAS CALCULATED IN * 

C * DEGREES FAHRENHEIT * 

C * SIGMAP IS THE VARIATION T N PHMAX * 

C * SIGMAS IS THE VARIATION ■ N CASE MATERIAL YIELD STRENGTH * 
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C * XI IS THE NUMBER OF STANDARO DEVIATIONS IN PHMAX TO BE USED * 

C * AS A BASIS FOR DESIGN ♦ 

C * X2 IS THE NUMBER ^F STANDARD DEVIATIONS IN SY TO BE USED AS ♦ 

C * A BASIS FO; JESIGN * 

C * SYCNOM IS THE NOMINAL YIELO STRENGTH OF THE CASE MATERIAL * 

C * IN LBS/INCH ♦ 

C * DCC IS THE ESTIMATED MEAN DIAMETER CF THE CASE IN INCHES * 

C * PSIC IS THE SAFETY FACTOR ON THE CASE THICKNESS * 

C * PELC IS TOE SPECIFIC WEIGHT OF THE CASE MATERIAL IN LBS/IN**3 * 

C . * LCC IS THE LENGTH OF THE CYLINDRICAL PORTION OF THE CASE * 

C * INCLUDING FORWARD AND AFT SEGMENTS IN INCHES * 

C * NSEG IS THE NUMBER OF CASE SEGMENTS * 

C * HCN IS THE AXIAL LENGTH OF THE NOZZLE CLOSURE IN INCHES * 

C * SYNNCM IS THE NOMINAL YIELO STRENGTH OF THE NOZZLE MATERIAL ♦ 

C * IN LBS/INCH * 

C * PSIS IS THE SAFETY FACTOR ON THE NOZZLE STRUCTURAL MATERIAL * 

C ' PS l A IS THE SAFETY FACTOR ON THE NOZZLE ABLATIVE MATERIAL * 

r * K1 AND K2 ARE EMPIRICAL CONSTANTS IN THE NOZZLE WT. EQUATION * 

** * PSIINS IS THE SAFETY FACTOR ON NOZZLE INSULATION * 

C * DEL INS IS THE SPECIFIC WEIGHT OF THE INSULATION IN LBS/IN**3 * 

;U01 CONTINUE 

C * KtH IS THE EROSION RATE OF INSULATION TAKEN CONSTANT * 

C ♦ EVERYWHERE EXCEPT AT THE NOZZLE CLOSURE IN IN/SEC * 

C ♦ KEN IS THE EROSION RATE OF INSULATION AT THE NOZZLE CLOSURE * 

C * IN IN/ SEC * 

C * OLINER IS THE SPECIFIC WEIGHT OF THE LINER IN LBS/IN**3 * 

C * TAUL IS THE THiCKNESS OF THE LINER IN INCHES * 

C * WA IS ANY ADDITIONAL WEIGHT NOT CQNSICEKED ELSEWHERE IN LBS * 


£ «4t***4t*4t4i* $***************'.:*¥¥***$***♦***** V***?*************#*-****** 

WRITE (6,610) DTEMPfSIGMAPtSI GMAS , X 1 , X2, SYCNOM, DCC ,PS I C, DELC, L 

ICC ,NSEG,HCN,SYNN0M,PSIS,P$IA,K1,K2,PS I INS, DEL INS,KfcH,KEN,L»L INER.TA 


2UL , WA 

832 CONTINUE 

► IF ( ISO.EQ.O .OR. Y.GT. 0.0) GO TO 80 

► RE AD ( 5,50b) PMOD,CMOD,PMU,CMU, ALPT S,T AUC 

*C * READ IN THE CONSTANTS FOR GRAIN DEFORMATION * 

► C ***NQTE*** NOT REQUIRED IF ISOO * 

► C * PM 00 IS THE ELASTIC MODULUS CF THE PROPELLANT IN PSIA I AT THE * 

►C * TEMPERATURE TORI * 

C * CMOD IS THE ELASTIC MODULUS CF THE CASE I ai THE BULK TEMPERATURE * 
‘ * OF THE GRAIN) IN LBS/IN**2 

►C * PMU IS POISSONS RATIO FOR THE PROPELLANT ( 0 ) THE BULK T EMPCKATURE* 
► r * OF THE GRAIN) * 

* CMU IS POISSON'S RATIO FOR THE CASE i .2 1 HE BULK TEMPERATURE OF * 
►C * THE GRAIN) * 

►C * ALPTS IS THE LINEAR COEFFICIENT OF THERMAL EXPANSION OF THE * 

►C * PROPELLANT IN IN/IN * 
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►C * TAUC IS THE CASE THICKNESS IN INCHES * 

►C a********************************** ****♦**♦••****•******♦***♦******♦* 

► WRITE (6, 799) PMOOt CMOD,PMU,CMU, ALPTS, TALC 

► 80 CONTINUE 

CALL AREAS 
IF(Y.LE.O.O) VC=VCI 
IFIABS(ZM).GT.O.O) GO TO 20 
*F ISUHAB.LE .0. 0) GO TO 31 
X=<ABPORT*-ABSLOTJ/SUMAB 

90 ’v%OZ=AT*X/APNOZ*t2.Ml.*BQT/2.*MNi*MNl)/TOP)**ZAP 
iF lABSIMNOZ-NNll.LE. 0*002) GO TO 2 
(L« l=MNOZ 
GO TO 90 

2 VNOZ=GAM*CSTAR*MNOZ*SQftT(U2./TOP)**lTOP/BOT)|/ll.+BOT/2.*MNUZ*MNO 
1Z) ) 

PRAT=I l» ♦BOT /2 •*MNOZ*HNOZ) *♦ (-GAM/BQT ) 

JROCK=AT /APNOZ 
SUMYA'-OELYM A6P2«-ABN2*ABS2 1 
IF(Y.EQ.O.O) SUHYA=0.0 
VC=VC+SUMYA 

► IF I ISO.EQ.O) GO TO 9 

► RATR=2.*//D0 ♦ OI/DO 

► RATR2=RATR**2 

► RATP=I2.*RATR2*< l.-PMU**2i/l 1.-RATR2)) /( PMGD*DO* ( l.-PMU*CMU)/ (CMOD 

► l*TAUC*2.l-PMU*l l.+PMU)*! 1 .-PMU**2 > * l 1 . *R ATK2 ) /I I. -RA TR2 ) ) 

► 9 CONTINUE 

IFIY.GT.O.O) GO TO 11 
€*♦** MOO I F ICAT ION MADE 1/7/76 

PONUZ=(QL*RHO«‘CSTAR*SUMAB/AT)**< U/ { 1. -Nil) *<1.«-(CAPGAM*JR0CK 1**2/ 
12. )**1N1/I 1.-N1)) 

IF C PONGZ.GT .P T RAN ) PONOZ 3 { Q2*RHU*CST AR* SUMAb/ AT )** 1 1 »/ I 1 .-N2 ) ) * 1 1 
l ( C APGAM* JROCK I **2/2. I **( N2/ ( 1.-N2)) 

MO I S=AT*PONOZ/CST AR 
P2=P0N0Z 
P0N0Z2-P0N0Z 
PNOZ=PRAT*PONOZ 

P4=2.*MDIS*VN0Z/l APHEAD+APNOZI+PNOZ. 

IF IGRA1N.E0.3) P4=MDIS * VNQZ/APNOZ ♦ PNCZ 
5 PNOZ= PRAT*PONOZ 

PHfcAD=2.*MOIS*VNOZ/( APHEAD+APNGZ ) +PNOZ 
IFIGRAIN.EU.3) PHEAO=MDIS * VNuZ/APNuZ ♦ PNOZ 
IF (PHEAD.LE.PTRAN|RHEA0=Q1*PHEAD**NI 
IFIPHEAO.GI . PTkAN) RHEAD=Q2 *PHEA0**N2 

zi t=mdi$*x/apnoz 

RN 1=RHE AO 
PHEA02=PH£AD 

C**»* MODIFICATION MADE 1/9/76 

3 IF ( ALPHA) 121,132,131 
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132 IF IPNOZ.LE.PTRAN) RN0Z=UI*PN0Z**N1 
IF(PNOZ.GT.PTRAN) RNQZ=Q2*PN0Z**N2 
GO TO 4 

131 IF (PNQZ.LE.PTRAN)RNOZ=RN1~((RN1-XG*PNOZ**XN-ALPHA*ZIT**.8/(L**.2*E 
IXP (BETA*RNl*RHO/ZIT) I )/( l.»ALPHA*ZI T**.d*bETA*RHO/ZlT/ (L**.2*EXP(B 
2ET A*RN1*RH0/Z I Til) ) 

IF IPNOZ.GT.PTRANIRNOZ=RNl-|(RN l-Q2*PNOZ*'*N2-ALPHA*Z IT**.8/CL**.2*E 
IXP(BETA+RNl*RHO/ZIT) ) )/( l.^ALPHA*ZIT**.8*BETA*RHQ/ZIT/(l.*'" Z EXPCB 
2ET A*RNl*RHQ/Z I T ) ) ) I 
IF (ABSIRN1-RNOZI.LE. 0.0021 GO TO 4 
RNi=RN0Z 
GO TO 3 

4 AVEi = (RHEAlHRNOZ 1/2. 

► IFIISO.EQ.1) GO TO 404 

► RH05=RH0 

► A8P0R6=ABP0RT 

► 404 CONTINUE 

IF(Y.GT.O.O) GO TO 7 

RN2=RNOZ 

RH2=RHFA0 

PONJ=PONOZ 

DPCDY=O.C 

AVE2=AVE 1 

7 RNAVE=(RNOZ<-RN2>/2. 

RHAV E= (R HEADER H2 )/ 2. 

► IF { I SO.EQ.O) GO TO 8 

► 401 CONTINUE 

► ETHETA=(PONOZ/PMOL>)*( RATP*C2.*(PMU**2-1. )/( l.— RATR2I ♦PMU*CMU*PMOO* 

► 100/(2. *CMQO* TAUC ) ) ♦( 1 . *R ATR2 ♦PMU* ( l .-KATK2-2 . *RATR2*PMU 11/ ( l.-RAlR 

► 22 ) ) 

► RH05=RH0/(( l. + ALPTSM TGR-TREF) )*( l.-PONOZM 1 .-2.*PMU )/ PMOD) )**3 

F IF (KOUNT.GE.l) GO TO 40b 

► ABP0R6-ABP0RT* (1.<-ETHETA*XETH) 

► 8 CONTINUE 

► IF ( PONOZ .L E. PT RAN )MGEN=RHOi>* CAVE 1*1 ABPOR 6+ABSLOT ) +Q1 ♦PONOZ **N 1*ABN 

k 10Z) 

IF (PGNOZ.GT.PrRAN)MGEN = RHOb*( AVE1M ABPCRt + ABSLCTI*U2*PONOZ**N2*ABN 
» 1UZ ) 

DRDY= ( A V£ 1— A VE 2 ) /DEL Y 
RdAR= (AVcl*AVE2)/2 . 

GM AX= 1 .0J02*MD 1 S 
GM IN- 0.999d*MD I S 
IF1Y.GT.0.0) GO TO 12 
GMAX= l.U01*MDI S 
GMIN=0.999*MDIS 

IF 1MGEN.GE.GMI N.ANO.MGEN.LE.GMAX) GC TO 6 
MO IS=MGEN 

PONQZ=MDIS*CSTAR/AT 
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IB 

16 

55 


C**** 


0ELY=ANS-YLED*0.05*DELTAY 

Y-ANS*0.05*DELTAY 


GO TO 5 

6 RE=2.*MDIS*X*L/IIAPNGZ*APHEAC»*MU> 

IF ( IGO.NE.O.AND. Y.LE.U.O) CALL IGNITN 
IF(Y.LE.O.Q) WRITE (6* 101 ) RE 
PONJ=PONOZ 
CALL OUTPUT 

10 IFIY.LE..05*TAUI GO TO 16 
SlNKl=VC/(CAPGAM*CSTAR|**2*RBAR*0PCDY/12. 

MASS=.01*M01S 
ANSA=Y*10.0*0ELTAY 
IF (KOUNT .GT.01 GO TO 16 

IF IABSISINKII.LE. MASS. AND. ANSA. Lt.ANS-XTJ GO TO 18 
GO TO 16 
0ELY=10.*DELTAY 
GO TO 55 
DELY=OEL TAY 
YLED= Y 
Y= Y+DELY 

AN S=T AU-A8S 11 / 2 ,) 

MODIFICATION MADE 1/9/76 
IF (Y.GE. ANS. AND. K HUNT .EQ.Q) 

IFIY.GE.ANS.ANO.KOUNT.EQ.OI 
0ELTAT=2.*DELY/IRHAVE*RNAVE» 

► DELTAT=OELTAT*( { l.-PONOZ*! l.-2.*PMU) /PMODI *1 1 . +ALPTS* ( TGR-TREF I i )* 

► l*3/{ i.+ETHETA) 

► DYETA=DELY*< II. +ALPTS* ( TGR-TREFI )*< l.-PCNGZ*! l.-2.*PMU) /PMGOl 1**3/ 

► lll. + ETHcTM 

». YE TA = YETA«-DYETA 

► IFIISG.EQ.OI YETA= Y 

SUM2- SUM A3 
RN2=RN0Z 

RH 2=RHEAD 
AVE2=AVE1 
GO TO 1 

11 MDIS=AT*P0N0Z/CSTAR 
GO TO 5 

12 DP CD Y = < PHEAD2 *P0NQZ2 I /( KNAVE *RHAVfc ) *DRDY ♦ l PH£ AD2 ♦PONUZ2 ) / 1 1 AUP2«-AB 
1N2 *ABS2 ) *2. J*DAUY 

IF (AHSIOPCDY).Gl.DPOUT.OR.Y.GE.XOoT I GC TC 25 

SIIMK1 = VC/(CAPGAM*C^TARJ**2*RGAR<-DPCDY/ 12 .♦ I PHEAU2*PGN0Z2 )/2.*(RNAV 
lE*KHAVE>/2. *lAbP2*ABN2*ABS2)/( 12.* I CS T AK *CAPGAM ) **2 ) 
STUFF=MG6N-SINKI 
MD I S= STUFF 
PONOZ-MDI S*CSTAR/AT 

IF ( Y.GE.0.9*( ANS-XT) )PUNOZ=Pl NJ*DPCDY* CC LY 
IF ISTUI F.GE .GMIN. AND. STUFF. LE.GMAX) GO TC 1A 
GO TO 5 
IA P1=P0N0Z 
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PONJ=PONCZ 

P0N0Z2=(Pl*P2)/2. 

P2=»f'0N0Z 

P3=PH6A0 

PHEAD2=(P3*P4l/2. 

P4=PHEA0 

HO IS=AT*P0N0Z/CSTAR 
OELTAT»2.*DELY/(RHAVE*RNAVE) 

» D£LTAT-DELTAT*(( l .-PONQZ* ( 1 .-2 .*PMU) / PMOO ) * ( 1. ♦ALPTS* ( TGR-TREF ) i i * 

► l*3/(i.*ETHETA) 

Z=Z*UELTATM RNAVE-RHAVE) 

T=T+DELTAT 

IF1Y.LT.ANS) CALL OUTPUT 
IFCY.LT.ANSI GO TO 10 

zw=z 

SUMBA= SUMAt 

P l=PONOZ 

RH2=RHEAD 

RN2=RNOZ 

RAVE=AVE1 

ABMAIN^SUMAB 

ABTG=0.0 

WR(TE(b.51) 

20 ANS2=TAU+ABS(ZW/2. ) 

KOUNT =KOUNT ♦ l 
IFIKOUNT.EQ.DCALL OUTPUT 
IF (KOUNT . E U. 1 ) GO TO 10 
OELYW=OELTAY 
0Y2=0ELYW 

IF(ZW) 32 * 32 » 33 

32 IFIY.LT.ANS2.AND.ABSIZWI.GT.DY2) GO TO 211 
SUMAB=ABMAIN+ABTT 

GO TO 31 

211 SUMDY=SUMOYfDELYW 

SUMAB=ll.+SUMDY/ZW-0ELYW/(2.*ZW) )*ABTO-l SUMDY/ZW-DELYW/(2.*ZW) )*AB 
IMA 1N+ABTT 
GO TO 31 

33 fF (Y.LT.ANS2.AND.ZW.GT.DY2) GO TO 21 
SUMAB=ABTO*ABTT 

GO TO 31 

21 SUMOY=SUMDY*DELYW 

SUMAB = ( 1 .-SUMD Y/ZW+DEL YW/ ( 2. *Z W) )*ABMA !N+( SUMOY/ZW-UELY*/ (2.*ZW) )* 
1A8T0+ABTT 

31 IF (SUMAB. LE. 0.0) PONUZ=PONOZ /2 . 

IF (SUMAB. LE. 0.0) GO TO 25 
MD I S=AT*PONUZ/CST AR 
AB AVE = ( SUMAB + SUMBA ) /2. 

SUMYA=DELY*ABAVE 
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VC=VC ♦SUMYA 

DADY= ( SUMAB-SUMBA) /OELY 
PBAR = (Pl+PONOZ 1/2. 

SUMBA=SUMAb 

22 IF ( PBAR.LE • P TRAN ID PC0Y=P6AR* CADY/ (1.— NU/ABAVE 
IF (PBAR„GT.PTRAN)DPCDY=PBAR*OAOY/( 1.-N2I/ABAVE 
PQN02=PONJ+OPCOY*OELY 

IF (PQNQZ .LE.O.U) PUNQZ=0.0 
IF (PONOZ.LE.PEXT) GO 10 25 
C **♦* MODIFICATION MADE 5/3/76 

IFIPONOZ.LE.PTRAN) Q2=Q1 
IF (PQNOZ.LE. PTRAN) N2=N1 
RNQZ=Q2*PGNOZ**N2 
C+*** MODIFICATION MADE 1/9/76 

RHEAD=RNOZ 

RBAR= (RHEAO+RAVEI/2. 

► IF ( I SO.EO. 0) GO TO 405 

► GU TO 401 

► 405 CONTINUE 

► M&EN=RHG5*(RN0Z*-RHEA0)/2.*SUMAB 

► IFUSO.EU.l) MGEN1*MGEN*I 1 . + ETHETA*XETH) 

► IF! ISO. EQ. 01 MGENl=MGEN 

GMAX= 1 .0002*M0 I S 
GMIN=Q.999d*MDlS 

SI NK^VC/ICAPGAM^CSTARI^^Z^RBAR^DPCDY/ 12 .+PBAR*ABAVE/ ( 12. *( CAP GAM* 
1CSTAR 1**2 )*RbAR 
STUFF=MGEN1-SINK1 
MOIS=STUFF 

1FISTUFF.GE.GMIN. AND. STUFF. LE.GMAX) GO TC 23 
PB AR = I P L +PONOZ 1/2. 

GO TO 22 

23 RHAVE=(kH2+RHEAC}/2. 

RNAVE = (Ki\2«-RNUZ1/2. 

RH2=RHEAD 

RN2=RN0Z 

PH £AD= PONOZ 

RAVE=RHF. AD 

P 1 = P()NOZ 

PQ NJ = PONOZ 

MDIS=AT*PONOZ/CSTAR 

IF ( ABS( DPCDY I . Gt.D p OUT) GO TO 25 

IF(Y.GE.XOUT) GO TO 25 

DELTA T = 2.*0f:LY/(KHAV£«-RNAVE) 

► DELTAT=0ELTAT*( ( U- PONOZ* I l.-2.*PMU)/PM0D)*I 1 . +ALPT5* l TGR-TREF 1 1 ) * 

» 1*3/1 l.fETHEJA) 

Z=Z+OELTAT*(UNAVE-RHAVEJ 

T=T*DELTAT 

CALL OUTPUT 
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GO TO 10 
25 RHEAD=0.0 
RN0Z=RHEA0 
PH£AD=P0N0Z 
MD IS=AT*PQNOZ/CSTAR 
WRITE (6*318) 

D£LTAT=2.*DELY/(RHAVE«-RNAVE) 

► OELTAT=OELTAT*( (l.-PONOZM l.-2.*PMU»/PMOD)*( 1 . *ALPTS* ( TGR-TREF ))) * 

► L»3/C l.+ETHETA) 

T* T+DELT AT 
CALL OUTPUT 
TIME=T 
0ELTAT=.5 
TIM=TIME*5. 

PHT=PHEAD 

SG=0.0 

29 T= T+OELTAT 

PHEAO=PHT/EXP( CAPG AM**2*AT*CSTAft/VC*( T-T IME) *12. ) 

PONOZ=PHEAD 
MD IS=PUi\UZ*AT/C STAR 
Y=Y+.5*RHEAD 
CALL OUTPUT 

IF(T.LT.TIM.AND.PHEAD.GE.0.04)G0 TO 29 
WP 1=G* SUMMT 
WP2=RHO*| VC-VCI )*G 
WP=( WP1+WP2 ) /2. 

► IF ( I SO.EQ. i I «P=WP1 
I SP= I TOT/ WP 

I SPVAC=l T VAC/ W P 
FAV= I TOT/T 
FVACAV=I TVAC/T 
PONAV=PONTOT/T 
LAM60 A= ( VC- VC I 1/VC 

WRITE (6, 1021 WP1 ,WP2,WP,PHMAX, ISP, ISP VAC, I IOT* IT VAC* F AV , FVACAV , PON 
1AV »VC ItVC»LAMBOA 
IF ( I wO . E 0 . 0 ) GG TO 903 

PMEOP=PHMAX*( 1 .♦Xl*bIGMAP ) *E XP ( P I PK*UT EMP ) 

SYC=SYCNQM*l 1.-X2*SIGMAS) 

TAUCC=PSIC*PME0P*0CC/I2.*SYC 1 

WCC=PI*T AUCC*DCC*0ELC*LCC*( l.+INSEG-l. 1 * ( AO. * TAUCC/LCC ) ) 
TAUCD=TAUCC/2. 

WC H-2. 5 *Pl/2.*0CC*0CC*TAUCD*DELC 

WC N=4.5*Pl/2.*0CC*HCN*TAUCD*DELC 

WC-WCC+wCH+WCN 

EPSIL = AE/PI/OT I/DT 1*4. 

WN=Ki*DTI*DT 1/ (1.+.5*SINIALFAN) )*( (EPSIL-SWRTIEPS1L ) l*PMEOP*DT l*PS 
11S/SYNN0M+K2*T*PSIA) 

WINS=T*PSI INS*DELINS*DCC*PI*lKEH*iDCC*.40*i S+NS) *TAU/2.*0. 15/ 
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IPSIINS*! LCC-TAU*! S*NS) ) ) *KEN*. 80*HCN) 

HL=T AUL*DL INER*Pl*DCC*( DCC/2 •♦LCC + HCN ) 

W ( ~WC+WN+H I NS+ WL+WA 
WM-W [ ♦ HP 
ZETAM=WP/WM 
RATIO 55 ITOT/WM 
WRITE!6,605) 

HR ITE ( 6* 601 ) PMEOP,TAUCC,wC,WN, WINS ,HL »W I ,WM , ZET AM, RAT 10 
903 CONTINUE 
NO UM= 1 

IF ( IPO.NE.O.AND.IPO.NE.2) CALL OUTPUT 
901 CONTINUE 

IF ( IOP.NE.G) CALL PLOT (0.0,0.0,999) 

STOP 

500 F0RMATI42X, 12) 

1903 FORMAT (6X,F6. 2, 10X.E11.4, 10X ,E 1 1. 4 , 8X, El 1.4 , / ,22X ,El 1 .4 ,9X, El 1 .4 ) 

11111 F0RMATI6X, I3.7X, 13 ) 

602 FORMAT ( 1H1 , 42X,21H.CCNFIGURAT ION NUMBER ,13) 

<,99 F0RMATI23F3.1) 

► 491 FORMAT ( 4X, l l ,9X, I 1,9X, 1 1 ) 

493 F0KMAT(4X,I1,15X,16I1,/,7X,I 1) 

► 492 FORMAT!//, 20X , 7H0PT l JNS, / , 13X. 5HIG0= ,11,/, 13X, 5HI W0= ,I1,/,13X,5H 

► *1S0= ,11) 

494 FORMAT (13X,5HIPO= , II,/, 13X, 12HNUMPLT ( JJ ) = , I 1 , 151 1H, , I 2 ) , 

2/, 13X , * I TEMP- • , 12 ) 

11112 FORMAT ! 13X, ‘NTA8= ' , l 3, / , 1 3 X , • NT ABY= ',13) 

► 505 FGRMAT(7X t E12.4,10X,C12.4,9X,F8.5»9X»Fb.5,/»10X,Ei2.4,lDX»F8.5) 

► 799 FORMAT!//, 15X.27HGRAIN DEFORMATION CON ST ANT S , / , 1 3X ,(>HP MOD 55 ,E15.4 

► 1 ,/ ,13X,6HCM00= ,E15.4,/, 13X,5HPMU= ,F8.4,/,13X,5HCMU= ,F8.4,/,13X, 

► 2 7H AL PTS= ,E15.4,/» 13X.6HT AUC 55 ,F8.4) 

501 FORMAT (7X,F10.0,/, 

2 4X,F9.6,3X,F7.5,3X,F6.3,6X,F5.2,5X,F6.2,4X,E11.4,/, 

26X ,F6.0) 

603 FORMAT! //, 20X , 26H PROPEL L ANT CHAR ACTER I S T ICS , / , 13X , 5HRH0= ,F9.6,/, 
113X,3HAl = ,F9.6,/,l3X,3HNl=,F6.3,/,13X, 7HALPHA= ,F6.2,/ , 13X,oHBETA= 
2 , F6.2,/,13X,3 MMU= , l PE l 1 . 4, / , 1 3X , 7HCSTAR= , 1PE 1 1 .*, /, 1 3X, • RN2N1 = 
2* , 1PE11.4) 

532 FORMAT (2X,Fd.2,5X,F6.2,4X,F7.2,5X,F6. 3,7X,F8.5,7X,Fd.5f/,l0X, 

1 F7.2,4X,F6.2,4X,Ft>.2,8X,F10./,6X,F8.2) 

604 FORMAT!// , 20X , 22 U3 AS I C MOTOR 0 I MENS I UNS , / , 1 3 X , 3HL= , F8 .2 , / , 1 3X , 5HT 
1AU= ,F6.2,/ , 13X, 4M0E = 

1 1 Pbll.4,/,13X, 5H0 T I - , l PE l 1 , 4, / , l 3X , 7HT FE T A= , 1PE1 1 .4 , / , 1 3 X , 7HALP 
3HA N= , 1 PE 11. 4,/, 13X,6HLTAP= , IPE 1 1. 4 , / , 13 X , 4HXT- , 1PL 1 1 . 4, / , 1 3 X ,4HZ 
40= ,1PE11.4,/,13X,BHCSTART- , 1 PE 1 1 .4 , / , 1 3X, 7HPTRAN= , 1PE1 1 . 4 , / , 1 3X 
5 , 4 FIN 2= , IPE 11. 4) 

10001 FORMAT (5X,F10.U) 

700 FORMAT (2F10.4) 

701 FORMAT 120X, • Y= • , IPEl 1 .4 , 10X , • TGR= *,IPE11.4) 
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503 FORMAT l7X,F6.3,5X,F7.2»7X,F7.2»7X,F5.4t3X,F6.2»3X,F8.0»/t 
15X,F7.4,8X,F8.5,5X,F8.2,7X,F7.3,5X,F7. 5, /, 5X , F 7. 3 , 5X,F 7.4, 
25X.F6.1) 

10002 FORMAT 1 13X» * TGR= *,1PE11.4) 

606 FORM AT(//»15X»27HB ASIC PERFORMANCE CONST ANTS, /, 13X ,8HDELTAY= ,F6.3 
l,/,i3X,6HXOUT= ,F8.2,/,13X,7HDP0UT= ,F 8. 1,/ , 13X, 7HZETAF= ,F7.4,/,1 

2 3X ,4HT B= » F6. 1,/,13X,4HHB= , F8.0,/, 13X ,5HGAM= ,F 7. 4, / , 13X , 7HERREF= 

3 ,F8.5,/,13X, 6HPRE F= ,F8.2 , / ♦ 1 3X , 7H0T R EF- ,F7.3 

4,/,13X,6HPIPK= ,F8.5,/, 13X,4HA2- ,F 8. 5 ,/ , 1 3X , 6HTREF = , F 7. 3» / 1 13X »6 
5HGAME= ,F7.4,/,13X,6HPEXT= ,F6.i) 

97 FORMAT! 3X,F 7. 1 , 5X, F6. 4.6X.F8. 1 , 7X,F 7. 1 , 7 X ,F7. I ,6X ,F5 .3 ,/ ,4X ,F5.2 t 
1 7X,F7.i,9X,F5.3, 7X,F7.3) 

842 FORMAT (20X.18HIGNIT ION C UNST ANT S ,/ , 13X ,4HKA= , F7. I « / , 13X, 4HK8= , 

1 F7.4,/, 13X,5HUf S= , Fd. L , / , 13X , 6HCS IG- ,F 7. 1 , / , l 3X , 6HPMIG= » 

2 F7. 1 ,/, 13X,5HTI 1= , F6. 3 , / , 13X , 5HT 12= , F5.2 ,/ , 13X, 6HRR IG= , 

3 F8.lt/tl3Xy 8HDEL T IG= , Fb. 3 , / , 1 3X ,6HPB I G= ,F7.3,//I 

600 FORMAT! 2 IX, F6. 2, 10X , F6. 3 , UX , F6 . 3 ,6X, F 5. 2, / , 5X , F 5. 2 , 10X.F 10 

I. 2 ,7X,F7. 2, 9X,F 5.2 ,8X,F6. 3,/ ,6X,F 8.2, 8X,F4.0, 7X,F /.2 , 10X,F10.2t8X, 
2F5.2t/t 7X.F5.2 ,6X,F7.4,6X,F7.4, 10X, F 5. 2. 10X, F 7.4, /, 6X, F 7.4, 7X ,F7.4 
3,10X»F7.4»8X,F7.4»6X,F9.2) 

610 FORMAT! 20X,19HINERT WEIGHT INPUTS,/ ti 3X, 

1 7H0TEMP* ,1PE11.4,/,13X,8HSIGMAP= , 1 PE 11 .4 ,/ , 13X, 8HSIGMAS= ,IPE11. 
24,/,13X,4HXl= .1PE11.4,/ ,13X,4HX2= , 1 PE 1 1 .4 , / , 1 3 X , 8HSYC.\iOM= ,1PE11 

3.4, /,13X,5HDCC= , 1PE 1 1 .4 , / , 1 3X , 6HPS I C= , IPE 1 1 . 4, / , 1 3X , 6HDELC= , IPE 

411.4, /,13X,5HLCC= , lPfc L l . 4 , / , 1 3X , 6HNSE G= , 1 PE 1 1 . 4 , / , 1 3X , 5HHCN= , IP 
5E1 1.4,/, 13X,8HSYNNUM= , 1 PE 1 1 . 4 , / , l 3X , 6hP S I S= , IPE 1 1 . 4 , / , 13X , 6HPS I A 
6= , 1PE11.4,/, 13X,4HK1= , 1 PE 1 1 . 4, / , 1 3X , 4HK2 = , IPE 1 1 . 4, / , 13X , 8HPS I I N 
7S= , 1PE1 1.4,/ , 13X, 8HUEL INS= , l PE 11 . 4 , / , 1 3 X ,5HKEH= , 1 PE 1 1 . 4, / , 1 3X, 5 
8HKEN= , IPE11.4,/.13X,8H0LINER= , IPE 11 . 4 , / , 13X , 6HT AUL= ,1PE11.4,/,1 
93X 4HWA= , IPE 1 1.4) 

101 FORMAT I// / , 33X , 29H ****$************************, / , 33X, 29H**** EQUI 
1LI BRI UM BURNING ***,/, 33 X , 29H** ***************************,//, 30X, 
225HINITIAL REYNOLDS NUM6£R= .1PE11.4) 

51 FORMAT137X, 2 3H**** **********«<‘**««‘***,/,37x,23H****TAIL OFF BEGINS 
1***#,/,37X,2 3H ****♦**♦«; *♦ fc*<t*#*#**$** f // j 
318 F0RMAT(37X,23H***********************,/,37X,23HBEG1N HALF SECOND T 
1 RA CE, /»37X, 2 3H ******* ************** **,//) 

102 FORMAT I13X, 5HWP 1 = , IPE 1 1.4 ,/ t i 3X,5H*P2 = , IPE 1 1 . 4 , / , 1 3X , 4HWP= ,IPEI 

II. 4, / ,13X,7HPHMAX= , IPE11.4,/,13X,5HISP= , 1 PE 1 1 . 4 , / , 1 3 X , 8HI SP VAC = 
2,lPEll.4,/,l3X, t.H I TO T = , 1PE1 1.4,/,13X, 7HITVAC= , l PE 1 1 . 4 , / , 1 3X , 5HF 
3AV= , 1PE11.4,/,13X,8HFVALAV= , IPE 1 1 .4 , / , 1 3X, 8HP0NAV= , iPEli.4,/,1 

3X,5HVCI= ,1PE11.4,/,13X,5HVCF= , IPEli.4,/,1 3 X , 8HL AMBDA-= .1PEJ1.4) 
605 FORMAT!///, 42X25HMGT0R WEIGHT CALCULATIONS) 

601 FORMAT! 13X,23HMAX EXPECTED PR£5SURE= , 1PE11.4,/ ,13X» 28HCYL 1 NOR I CAL 
1 CASE THICKNESS^ , IPE 1 1 . 4 , / , 1 3X , 9HC A SE nl = , 1 PE 1 1 . 4 , / , 1 3X , 1 1HN0ZZL 
2 E rtT = ,1PE11.4,/,13X, 15HINSULATICN WT= , IPE 1 1 . *, / , 1 3X, 10HL I NER WT= 
3 , lPEli.4,/,13X,16HTUTAL INERT WT= , 1 PE 1 i .4 , / , 13X , 2 OHT OT AL MOTOR W 
4E l GHT = , LPE11.4,/, 13X,7HZETAM= , IPE 1 1 .4, /, 13X , 21 HRAT 10 OF I TOT TO 
5WM = , IPE 11.4) 

END 


- 94 - 



Table C-3 (Cont'd) 


SUBROUTINE AREAS 

C **************************** ***************************************** 

C * SUBROUTINE AREAS CALCULATES BURNING AREAS AND PORT AREAS FOR * 

C * CIRCULAR PERFORATED IC.P.) GRAINS AND STAR GRAINS OR FOR A * 

C * COMBINATION OF C.P. AND STAR GRAINS * 

C ********************************************************************* 

INTEGER STAR, GRAIN, ORDER, COP 

REAL MGEN,MDIS,MNOZ,MNi, JROCK , N, L , ME 1 , ME , ISP, I TOT, MU, MASS, ISP VAC 

REAL LGC I , LGN I , NS , NN, NP, LGSI »NT»LTP»LGC»LS,LF 

REAL M2 , MDBAR, I SP2 , I TVAC , L I • L 2 , LFW.LFW SQO 

COMMON/C ONST1 /ZW,Afc, AT, THETA, ALFAN 

COMMON/CONST 3/ S, NS, GRAIN, NT AtlY.NCARO 

COMMON/C 0NST4/DELD I , DO,ZO,DI 

COMMON/ VAR IAL/Y»T»DELY»DELTAT, PONOZ »PHEAO »RNOZ ,RHEAD» SUMAB, PHMAX 
COMMON/ VAR I A2/ ABPORT , AdSLOT , ABNOZ , APHE AO , APNOZ ,0 ADY , ABP2 , AbN2 , AbS2 
COMMON/VARI A3/I TOT , IT VAC , JROCK, I SP, ISPVAC,MDIS,MNUZ, SG,SUMMT 
COMMON/ V AR I A4/ RNT ,RHT»SUM2»RI,R2»R3,RHAVE*RNAVE»RBAR*YB*KQUNT ,TL 
COMMON/VARI AS/ ABMA I N , ABTO ,SUMOY , VC I , AB TT , PTR AN 
COMMON/VAR I A8/YDI 

t COMMON/VARI 20/YETAtXETH, ISO 

DATA PI/3.14159/ 

ABPC=0.0 
ABNC=0.0 
ABSC=0.0 
AB PS=0 .0 
ABNS=0.0 
AB SS=0.0 
DA dT=0.0 
SG=0.0 
VC IT=0.0 
ANUM=PI/4. 

P I D2=P 1/2. 

RNT=RNT*RNOZ*DELTAT 
RHT=RHT+kHEAU*DELT AT 
IF(Y.LE.O.O) AGS=0.0 
K= 0 

IF(ABSIZW).GT.O.O) K=l 
YB=Y 

IF ( K. EQ. L ) Y= Y B-SUMDY/2 . 

2 IFIK.EQ.2) Y=YB+ABS(Zri)/2.-SUMUY/2. 
v IF I ISO.EQ.O) YE T A= Y 

r r (Y.LE.O.O) READ I 5, 500) I NPUT , GRA I N , ST AR , NT , CRDE R , COP 

C ********************************************************************* 


C * RfcAU THE TYPE OF INPUT FOR THE PROGRAM AND THE BASIC GRAIN * 
C * CONFIGURATION AND ARRANGEMENT * 
C * VALUES FOR INPUT ARE * 
C * 1 FOR ONLY TABULAR INPUT * 
C * 2 FOR ONLY EQUATION INPUTS (EQUATIONS ARE BUILT * 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


* INTO THE SUBROUTINE I 

* 3 FOR A COMBINATION OF 1 ANC 2 

* VALUES FOR GRAIN ARE 

* l FOR STRAIGHT C.P. GRAIN 

* 2 FOR STRAIGHT STAR GRAIN 

* 3 FOR COMBINATION OF C.P. AND STAR GRAINS 

* VALUES FOR STAR ARE IWAGON WHEEL IS CONSIDERED A TYPE OF 

* STAR GRAIN IN THIS PROGRAM) 

* 0 FOR STRAIGHT C.P. GRAIN 

* 1 FOR STANDARD STAR 

* 2 FOR TRUNCATED STAR 

* 3 FOR WAGON WHEEL 

* VALUES FOR NT ARE 

* 0 IF THERE ARE NO TERMINATION PORTS 

* X WHERE X IS THE NUMBER OF TERMINATION PORTS 


* 

♦ 

* 

* 

* 

* 

♦ 


VALUES OF ORDER ESTABLISH HOW A COMBINATION C.P. AND STAR 
GRAIN IS ARRANGED 

1 IF DESIGN IS STAR AT HEAO END AND C.P. AT NOZZLE 

2 IF DESIGN IS C.P. AT HEAD END AND C.P. AT NOZZLE 

3 IF DESIGN IS C.P. AT HEAD END AND STAR AT NOZZLE 

A IF DESIGN IS STAR AT HEAD END AND STAR AT NOZZLE 

♦♦•NOTE*** IF GRA IN=1 « VALUE CF ORDER MUST BE 2 


LOOO CONTINUE 


* 

* 

* 

* 

* 

* 

* 

* 


♦♦♦NOTE*** IF GRAIN=2» VALUE CF ORDER MUST BE 4 
VALUES FOR COP ARE (APPLICABLE TO C.P. GRAINS ONLY) 

0 IF BOTH ENDS ARE CONICAL CR FLAT 

1 IF HEAD END IS CONICAL OR FLAT AND AFT END IS 

HEMISPHERICAL 

2 IF BOTH ENDS ARE HEMISPHERICAL 

3 IF HEAD END IS HEMISPHERICAL ANO AFT END IS 

CONICAL OR FLAT 


* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


it******************************************************************** 


IF(Y.LE.O.O) WRITE (6»607 ) 

IF(Y.LE.O.O) WRI T£ ( 6 » 600 ) I NPUT , GRA l N , £T AR , NT , ORDER , COP 
IF ( INPUT .EQ.2 ) GO TO 12 
IFIY.LE.O.O) GO TO 6 
IF (K. EQ.2) GO TO 91 
IF (K.EQ. 1)Y=YB 
IF(YT.LE.Y) GO TO 8 
9 DE NOM= YT-YT2 

SLOPE 1= l ABPK— ABPK2 ) /DENUM 
SLOPE 2= ( ABSK— ABSK2 ) /DENOM 
SLOPE 3= ( ABNK- ABNK2 ) /DENOM 
SLOP £4= ( APHK— A PHK2 ) /DENOM 
SL0PE5= ( APNK-APNK2 )/ DENOM 
B1=ABPK-SL0PE1*YT 
B2=ABSK-SLOPE2*YT 
B3=ABNK-SLOPE3*YT 
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B4=APHK-SL0PE4*YT 

B5=APNK-SL0PE5*YT 

ABPT=SLOPEi*Y*Bl 

ABST=SLOPE2*Y+B2 

ABNT=SLOPE3*Y+83 

APHT=SlOPE4*Y+B4 

APNT=SL0PES*Y«-B5 

YB=Y 

IFIK.EQ.l) Y=YB-SUMDY/2. 

91 IF ( INPUT. EQ. 3) GO TO 3 • 

GO TO 52 

6 RE ADI 5» 507) YT , ABPK , ABSK, ABNK, APHK , APNK , VC I T 

NC ARD=NC ARD+- 1 

C ************************♦****************♦*♦******¥*********** ******* 


C * READ IN TABULAR VALUES FUR Y=0.0 (NOT RtUUIRED IF INPUT=2) * 
C * * 
C * ABPK IS THE BURNING AREA IN THE PORT IN IN**2 * 
C * ABSK IS THE BURNING AREA IN THE SLOTS IN IN**2 * 
C * ABNK IS THE BURNING AREA IN THE NOZZLF END IN IN**2 * 
C * APHK IS THE PURT AREA AT THE HEAD END IN IN#*2 * 
C * APNK IS THE PORT AREA AT THE NOZZLE END IN IN**2 * 
C * VC I T IS THE INITIAL VOLUME UF CHAMBER GASES ASSOCIATED WITH * 
C * TABULAR INPUT IN IN**3 * 


C *** ****************** ***** ********************* *************** ******* 
WR I TE I 6 » 6 1 0 I 

WRITE ( 6 i 583) ABPK , ABSK , ABNK , APHK , APNK 
WR I T£ ( 6, 584 ) VC l T 
AB PT= ABPK 
ABST=ABSK 
ABNT = ABNK 
AP HT = APHK 
APNT=APNK 
YT 2= Y T 

IF ( INPUT. EQ. 3) GO TO 3 
VC 1 = VC I T 
GO TO 52 
8 Y T ? = Y 1 

ABPk2=ABPK 
ABNK? = ABNK 
ABSK? = ABSK 
APHK2= APHK 
> ? NK 2 = APNK 

READ l 5 » 5 Q5 ) YT , ABP K , ABSK , ABNK , APHK , APNK 
NC ARU = NCARD«-l 

**<t*********** V 3*******V**<>******************^*********************** 

* READ IN TABULAR VALUES FOR Y=Y (NOT REQUIRED FOR INPU T=Zi * 

* (NOTE THAT TABULAR VALUE CARDS FOR Y GT 0 CO NOT IMMEDIATELY * 

* FOLLOW THOSE FOR Y EQ 0 IN THE DATA DECK! * 
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C **************** **************^****** **************** ************ 

WRITEl6,6ll) YT 

WRITE (6, 583 J ABPK. ABSK, ABNK , APHK , APNK 
GO TO 9 
12 ABPT=0.0 
A8NT=0.0 
A8ST=0.0 

3 IFIGRAIN.NE.2) GO TO 4 
ABPC=0.0 

•>ABNC=0 .0 
A8SC=0.0 
GO TO 7 

4 IF (Y.LE.O.O) READ! 5*5011 DO, 01 , OELD I , S ,THETAG, LGC l , LGN I , THE TCN, THE 
ITCH 

C ^^^*^^irn^t****************** ****** *v**** ***************************** 

C * RE AO IN BASIC GEOMETRY FOR C.P. GRAIN INGT REQUIRED FOR * 

C * STRAIGHT STAR GRAIN] * 

C * DO IS THE AVERAGE OUTSIDE INITIAL GRAIN DIAMETER IN INCHES * 

C * 01 IS THE AVERAGE INITIAL INTERNAL GRAIN DIAMETER IN INCHES * 

C * OELOl IS THE DIFFERENCE BETWEEN THE INITIAL INTERNAL GRAIN * 

C * DIAMETER AT THE NOZZLE END OF LGCI AND DI IN INCHES * 

C * S IS THE NUMBER OF FLAT BURNING SLOT SIDES (NOT INCLUDING * 

C * THE NOZZLE END) * 

C * THETAG IS THE ANGLE THE NOZZLE ENO CF THE GRAIN MAKES WITH * 

C * THE MOTOR AXIS IN DEGREES * 

C * LGCI IS THE INITIAL TOTAL LENGTH OF THE CIRCULAR PERFURAT ION * 

C * IN INCHES * 

C * LGNI IS THE INITIAL SLANT LENGTH OF THE BURNING CONICAL * 

C * GRAIN AT THE NOZZLE END IN INCHES * 

C * THETCN IS THE CONTRACTION ANGLE OF THE BONDED GRAIN IN DEG. * 

C * THETCH IS THE CONTRACTION ANGLE AT THE HEAU END IN DEGREES * 

C ****<t#***##**$<:>Mt*$*<t**#********#**$i|r*<:***$*<i<:*#**j{<*$*<M(t<c#*$**#<s*#*<:$ 

IF { Y.LE.O.O) WRITE (6,601 ) DO , D I , UELD I , S , THET AG , LGC I , LGNI , THETCN, TH 
1ETCH 

IF (Y.LE.O.O) THETAG=THETAG/57.2957B 

IF ( Y.LE.O.O) THE TCN= THETCN/ 5 7. 2957 8 

IF ( Y.LE.O.O )THETCH=ThE TCH/ 5 7. 2957b 

D0SQD=00*U0 

DI SOD=DI*Ol 

BNUM=ANUM*DOSQD 

TLL^TL 

<F { ORDER. GE. 3) TLL=0.0 

YDI=2.*Y+DI 

YD I SQD= YD I #YD I 

AB SC= S*A NUM* ( OOSQU— Y D I SOD ) 

IF (ABSC.LE.O.O) ABSC=0.0 
IF (YDl .GT.DO) GO TO 100 
IF1THETAG.GT. 0.08727) GO TO 101 
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IFICOP.EQ.O) GO TO 700 
IFCCQP.EQ.l) GC TO 701 
IFJC0P.EQ.2) GO TO 702 
CHCK l-DQSQD— YD l SQO 
IF (CHCKl «L T .0. 0 ) CHCK 1=0. 0 

LGC=LGCI-C SORT ( DOSQD-DISQD J-SQRT ( CHCK l ) I /2.- Y*COTANC THETCNI 
GO TO 710 

702 CHCKl=DOSQD-YDISQD 

IFICHCKI.LT. 0.0) CHCK1=0.0 
IFCCHCKi.LT. 0.0) CHCK1=0.0 
LGC=LGC 1 -( SORT ( DOSQD-DISQD )-SQRT CCHCK l ) ) 

GO TO 710 

701 CHCK2=0QSQD-(Y0I*DELD!)**2 
IF CCHCK2.LT. 0.0) CHCK2=0.0 

LGC=LGC I - ( SORT ( DOSQO- lDI‘t-DELDI)**2) — SQRT (CHCK2) )/2. 
i-Y*COTANITHETCH) 

GO TO 710 

700 LGC=LGCl-Y*(COTANl THETCN)+COTAN(THETCH) ) 

► 710 ABPC=Pl*YDi*lLGC-TLL-S*YETA) 

A8NC=0.0 
GO TO 732 
101 CONTINUE 

IFICOP.EQ.O. OR. COP. EQ.i) GO TO 720 
CHCK 1=00 SQD- YD I SQD 
IFCCHCKI.LT. 0.0) CHCK1=0.0 

A8PC=PI*Y0I*( LGCI- ( SORTC DOSQD-Dl SQD ) -SQR Ti CHCK 1) )/2.-TLL 

► 2-( S+TAN( THETAG/2.) J*YETA) 

GO TO 730 

► 720 A8PC=PI*YDI*(LGCI-Y*C0TANl THETCHJ-TLL- l S+TANC THETAG/2. ) ) *YET A) 

730 IF (COP.EQ. l.OR. COP . EQ. 2) GO TO 731 

A8NC=PI* (LGNI -Y*CO TAN l THE T AG + T HETCNJ-Y * TAN ( THETAG/2. ) )*( DI + 

1 DELDH-Y*LGN1*S1N l T HET AG ) ♦ Y* S l N l THETCN) /S IN ( THE TAG+THE TCN) ) 

GO TO 732 

731 IF (Y.LE.O.O) GO TO 7311 
GO TO 7312 

7311 R7 = ( (DI+DELOl ) /2 . *1 GN I * S I N ( T HET AG ) >*CQSl THETAO-SINC THETAG) * 

1 iQRT((DU/2.)**2-( ( 0 I ♦ 0 F L D 1 ) / 2 . + LGN I * S I N { THE T AG ) )**2) 

7312 IF (R7 + Y.LT. I 00/2. ) *COSl T HE TAG) ) GO TO 1 i 111 
A8NC=Pl*ILGNl «• ( l./SINC THE TAG) )*( (DU/2. )-LGNI *S INC THE TAG) 

l- 1 GI+DELDI )/2. )-Y*CO TAN ( T HE 1 AG )- Y* T AN ( THETAG/2.) )*( IDI+DELDI ) 
2 / 2 .* Y*D0/2. ) 

GO TO 22222 

lilll RPR=SQRT( ( (D0/2.)**2)~R7**2)-SQRT1 I ( D0/2.)**2)-(R7*-Y)**2) 

A3NC=PI*(LGNI-KPR— Y*TANl THE T AG/2 .))* ( 1 1)1 +DELDD/2 .+SQRTC (DO/ 
l 2. ) **2- (R7 + Y ) **2 ) *SIN( THET AG ) + Y +1 R 7+ Y ) *CUS( THE TAG) ) 

22222 CONTINUE 

732 IF (A8PC.LE.0.0) ABPC = 0.0 
IF (ABNC.LE .0.0) ABNC^O.O 
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GO TO 5 
100 ABNC-0.0 
ABPC=0.0 
5 DH=DI-ZO 

APHT=ANUM* ( 0H»2.»RHT ) **2 
IF CAPHT.GE.BNUM) APHT=BNUM 
IFCK.LT.2) APHT l=APHT 
APNT=ANUM* (01 *DELD I f 2«*RNT ) **2 
IFIAPNT .GE.BNUM) APNT=BNUM 
IFCGRAIN.NE.il GO TO 7 
ABPS=0.0 
ABSS-0.0 
ABNS=0.0 
GO TO 50 

7 IF ( Y.LE.O.O) RE AD( 5 » 502 I NS , LGSI , NP , RC, F ILL, NN 
***** ********************** ************** ******************* ******** 

* READ IN BASIC GEOMETRY FOR STAR GRAIN C NOT REQUIRED FOR 

* STRAIGHT C.P. GRAIN) 

* NS IS THE NUMBER OF FLAT BURNING SLOT SIOES (NOT INCLUDING 

* THE NOZZLE END) 

* LGSI IS THE INITIAL TOTAL LENGTH OF THE STAR SHAPED 

* PERFORATED GRAIN IN INCHES 

* NP IS THE NUMBER OF STAR POINTS 

* RC IS THE AVERAGE STAR GRAIN OUTSIDE RACIUS IN INCHES 

* FILL IS THE FILLET RADIUS IN INCHES 

* NN IS THE NUMBER OF STAR NOZZLE END BURNING SURFACES 
********** ******************* ********** **** **** ********************* 

If(Y.LE.O.O) HR I TE (6.602 ) NS, LGSI ,NP ,RC , F I LL ,NN 
PI DNP=PI /NP 
RC SQO=RC*RC 

► IFCISO.EQ.l) Y E=Y 

► IFCISO.EQ.l) Y=YETA 
FY=FILL*Y 

FY SQO=FY*F Y 

IFCSTAK.EQ.il GO TO 20 
IFCSTAR.6Q.2 ) GO TO 201 
IF fY.GT.O.O) GC TO 179 

ADC « ,421) TAUWW.Ll ,L2,Al PHA 1 , ALPHA 2 » HH 
C *********************v*»v ******************************** *********** 
C * READ IN GEOMETRY FOR WAGON WHEEL C NCT REQUIREO FOR STANDARD 
Z * OR TRUNCATED STAR GRAINS) 

* TAUWW IS THE THICKNESS OF THE PROPELLANT WEB IN INCHES 

C * LI AND L2 ARE THE LENGTHS OF THE TWO PARALLEL SIOES OF THE 

: * TWO SET " OF STAR POINTS IN INCHES 

C * AL PHA 1 AND ALPHA2 ARE THE ANGLES BETWEEN THE SLANT SIDES OF 

C * THE STAR POINTS CORRESPONDING TC LI AND L2, RESPECTIVELY, 

C * AND THE CENTER LINES OF THE POINTS IN CEGREES 

T * HH IS HALF THE WIDTH OF THE STAR POINTS IN INCHES 
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C ********************************************************************* 
WRITE (6.422 ) T AUWW »L l *L2 1 ALPHA l » ALPHA2 »HW 
ALPHAl=ALPHAl/57. 29578 
ALPHA2*ALPHA2/ 57.29578 
ALP2*ALPHA2 
XL 2-L 2 

LF W=RC— T AUWW—F ! LL 
LF WSQO=LF W*LFW 
THETFW=ARSIN(( HW*F ILL )/LFW) 

SLFW=LFW*SINlTHfcTFW) 

179 KKK=0 
SG-0.0 

ENUM= ( RCSQO-LFwSQD— FYSQD )/ I 2 • *LF W*F Y ) 

AL PHA2-ALP2 
L2=XL2 

190 YTAN=Y*TAN( ALPHA2/2.I 
COSALP=COS( ALPHA2) 

S I NAL P=S I N ( AL PHA2 ) 

IFIYTAN.GT.L2) GO TO 182 
IF (FY.GT.SLFW) GO TO 181 

SGW-NP*! L2— 2.*YTAN*( SLFW-F I LL )/S INALP- Y*COT AN ( ALPHA2 ) *FY* 

1 IPID2+THETFW)+(LFW+FY)*‘C PI ONP—THETF W ) ) 

GO TO 183 

181 IFIY.GT.T AUWW ) GO TO 184 
SGW-NP*lFY*(PlDNP»ARSlN(SLFW/FY))*»2IDNP-THEr« : W)*LFW) 

GO TO 183 

184 $G W=NP*F Y* ( THE TFW+ARS INI SLFW/FY )-ARCOS I E MJM) ) 

GO TO 183 

182 YPO=— SLF W 

IFIALPHA2.GE.PID2) GO TO 222 
Q=-FILL+L2*TAN(ALPHA2)-Y/CCSALP 

XP I = ( — Q*T AN ( AL PHA2 )-SQRT (-Q*U+FYSQD/CCSALP*COSALP ) )*COSALP*COSALP 
YP I = XPI*TAN( ALPHA 2 ) *-Q 
XPO=< YPO-Q) *COTAN( ALPHA 2) 

GO TO 223 

222 XP 1 = Y-L2 

YP I SUR T ( F YSQD-XP I*XPI ) 

XPU=XPI 

223 FYLS = SQRT(SLFH*SIFW + XPI*XPU 
XPI02=I XPi-XPOI*( XPl-XPO) 

YP 102= IYPI-YPG)*! YPI-YPO) 

»F (FY.GT.FYLS) GO TO 186 
IF(Y.GE.TAUWW) GO TO 185 

SGW=NP*(SURT (XP102+YPI02 ) *F Y *1 P I D2+THE TF W-AR S IN (XP I /F Y ) )+(LFW*FY)* 
1 IPIONP-THETFW) ) 

GO TO 183 

185 SG W=NP*t SQR T(XP102*YPIU2)*FY+CPI D2-ARS INIXPI/FY )— ARCOS I ENUM ) ) ) 

GO TO 183 
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186 IFIY.GT.TAUWW) GO TO 187 
SGW=NP*(FY*(PIDNP*ARSlNISLFw/FYI M-I PI DNP-THE TFwl *LFW ) 

GO TO 183 

187 $GW=NP*FY*(THETFWfrARS!N(SlFW/FYl-ARCOSlENUMl 1 
183 IF J SGW.LE.0.0 ) SGW-0.0 

IFIY.GT.G.Oi GO TO 188 

AGS2* • 5*1 PI*RC SQD-NP*LF W*SlFW* ICOSCTHETFW)-SIN(THETFW) *COT AN ( ALPHA 
L 2)-2.*(L2*FlLL*TANI ALPHA2/2.) J/LFWl-I Pl-THETFW*NP|*LF WSQD-2. *NP*F 

2 lLL*(L2*SLFW/SINAL0*LFw*(PlONP-THETFWl+CPIDNP+PID2-l./SINALP)* 

3 FILL/2. U 
AGS s AGS+AGS2 

188 CONTINUE 
SG=SG*SGW 

IFIKKK.EQ.il GO TO 24 
L2=L1 

AL PHA2=ALPHAi 
KK K= 1 
GO TO 190 

201 IF(Y.LE.O.O) READ! 5»503) RP,TAUS 

C ********************************************************************* 


C ♦ READ IN GEOMETRY FOR TRUNCATED STAR (NCT REQUIRED FOR ♦ 
C * STANOARD STAR OR WAGGN WHEEL ) * 
C * RP IS THE INITIAL RADIUS CF THE TRUNCATION IN INCHES * 
C * TAUS IS THE THICKNESS OF THE PROPELLANT WEB AT THE BOTTOM * 
C * OF THE SLOTS IN INCHES ♦ 


C ******** ************************************************* ************ 

IF(Y.LE.O.O) WRITE (6*603 ) RP,TAUS 

THETAS=PI DNP 

RPY=RP+Y 

LS=RC-TAUS-F ILL-RP 
RPL=RP*LS 

THETS1=THETAS-ARSIN(FY/RPY1 
IF (THETS1.L6.0.0) GO TO 110 
IF I Y. LE. TAUS 1 GO TO 103 

THET AC-ARS IN ( I RCSQD-RPL*RPL-FY SQD ) / I 2 . *F Y*RPL ) ) 

IF (THETAC.GE.O.O) GO TO 104 
IF (Y.LT.RC-RPI GO TO 105 
SG =0 . 0 
GO TO 14 

103 SG=2.*NP*(RPY<‘THETS1*LS-IRPY*CQSI THETAS-THETSl I-KP ) *P I U2*F Y ) 

GO TO 14 

104 SG=2.*NP*(RPY*THETS1*LS-(RPY*C0S(THETAS-THETSII-RP)*FY*THETAC) 

GO TO 14 

105 SG=2.*NP*UPY*THETSi + SURT(RCSQ0-FYSQD)-Si.kT< RPY+RPY-FYSQD) ) 

14 IF (Y.L6.0.0) AGS=P I * I RE SQD-KP*KP |-NP* (PI^FILL^F ILL/2. ♦2.*IS*FILL) 
GO TO 31 

110 THETAF=THET AS 

THETAP=2.*THETAS 
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TAUWS=TAUS 
GO TO 111 

20 IF IY.GT.0.0) GO TO 1791 

READ! 5 • 504 I THE TAF , THETAP, TAUWS 

C ********************************************************************* 


C * READ IN GEOMETRY FOR STANOARD STAR I NOT REQUIRED FOR * 
C * TRUNCATED STAR OR WAGON WHEEL) * 
C * THETAF IS THE ANGLE LOCATION OF THE FILLET CENTER IN DEGREES * 
C * THETAP IS THE ANGLE OF THE STAR POINT IN DEGREES * 
C * TAUWS IS THE WEB THICKNESS OF THE GRAIN IN INCHES * 


C ******************************* +** *********************** ************ 

WRITE ( 6 * 604 ) THE TAF , THE TAP , TAUWS 
THETAF=THETAF/ 57.29578 
THETAP=TH£TAP/57.29578 
THETAS=P!/NP 
THETSl=l.OO 

111 lf=rc-tauws-fill 

1791 CNUM=(Y*FILU/LF 

DNUM= S I N ( THE TAF ) /S INI TH6TAP/2. ) 

ENUM= ( RCSQU-LF *LF-F YSQD l/I2.*LF*FY) 

FNUM=S INI THETAF )/CCS( THE TAP/ 2. ) 

IF ICNUM.LE.FNUMJ GO TO 106 
IFlY.LE.TAUwS) GO TO 107 

SG=2.’*NP*FY*t TH£TAF*ARSIN|SINITHETAF)/CNUM)-ARCUSIENUM) ) 

GO TO 23 

106 IFIY.LE. TAUWS) SG=2. *NP*LF *1 DNUM + CNUM* I P IC2* THETAS-THET AP/2 . 
1-COTANI THE TAP/ 2.) )*T HE TAS- THETAF) 

IFIY.LE. TAU»S> GO TO 23 

SG=2. *NP*I FYMARSlNl ENUM ) *T HE TAF— THE! AP/2. )*LF*ONUM-FY*COTANI THETA 
IP/2.) ) 

GO TO 23 

107 SG=2.*NP*LF*ICNUM* I THETAS* ARSi N I S IN l THETAF ) /CNUM) ) *THE TAS-T hETAFI 

23 IF ITHETSl.LE.O.O) GO TO 1* 

IFIY.LE. 0.0) AGS=PI*RC**2-NP*LF*LF*tSIN( THETAF I* I COS I THETAF )- 
1SI NITHETAF>*COIAN( THETAP/ 2.) ) *THETAS-THE TAF*2.*F ILL/LF*! SIN I THETAF 
2 I /SIN l THE T AP/2 .) ♦ T HE T AS- T PET AF *F ILL/I2.*LF ) * I P I D2*THETAS-THE 
3TAP/2.-C0TAN ( THETAP/2. ) ) )) 

24 CONTINUE 

31 IF ISG.LE .O.U) SG =0.0 

IF (K.EQ.0.0R.K.EQ.2) SGN=SG 

IFIK.LE.l) SGH=SG 

IFIY.LE. 0.0) SG2=SG 

IFIK.E0.2) GO TO 37 

RAVEDT=R 1*1 SG*SG2)/2.*K6AR*DELTAT 

RN0T=R2*ISG*SG2)/2. * RNAVE*OE L T AT 

RHQT=R3*t SG*SG2)/2.*RHAVE#DELTAT 

R1 =RAVEDT 

R2=RNDT 
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R3=RHDT 
SG2=SG 
GO TO 38 

37 IF IKOUNT .NE. 1 ) GO TO 39 
SG 3-SG 

R4=Rl 

R5=R2 

R6=R3 

39 RAVE0T=R4*ISG+SG3)/2.*RBAR*DELTAT 
RNDT=R5*(SG*SG3)/2.*RNAVE*DELTAT 
RHDT=R6M SG+SG3 )/2«*RHAVE*OELTAT 
R4=RAVEDT 
R5=RNDT 
R6=RH0T 
SG3= SG 

38 AB SS= ( AGS-RAVEDT ) *NS 

IF (ABSS.LE.O.O.OR. SG.LE.O.O) ABSS=0.0 
ABNS= I AGS-RNDT )*NN 

IF (ABNS.LE. 0.0. OR. SG.LE.O.O) ABNS=0.0 
IF l OROER .LE. 2 ) ABPS= I LGS I-Y* ( NS*NN) ) *SG 
IF IOROER.LE.2 ) GO TO 36 
A8PS= I LGSI-TL-Y*! NS*NN) )*SG 
36 PI RCRC^P I*RCSQD 

APHS=P I RCRC- AGS+RHDT 

IF l AP HS.GE.PIRCRC.GR. SG.LE.O.O) APHS-PIRCRC 
APNS-P IRCRC-AGS+RNDT 
IF IK.LT.2) APHS1- APHS 
IF(APNS.GE.PIRCRC) APNS=PIRCRC 
► IF I I SO.EQ. 1 ) Y=YE 

50 IF(NT.EQ.O.O) GO TO 371 

IF (Y.LE.O.O) REAOI 5» 506) L TP , OTP , THETTP , TAUEF F 
C ******************************* ************************************** 


C * READ IN GEOMETRY ASSOCIATED WITH TERMINATION PORTS (NOT * 
C * REQUIRED IF NT = 0) ♦ 
C * LTP IS THE INITIAL LENGTH OF THE TERMINATION PASSAGES * 
C * IN INCHES ♦ 
C * DTP IS THE INITIAL DIAMETER OF THE TERMINATION PASSAGE * 
C * IN INCHES * 
C * THETTP IS THE ACUTE ANGLE BETWEEN THE AXIS OF THE PASSAGE * 
r * AND THE MOTOR AXIS IN DEGREES ♦ 
C * TAUEFF IS THE ESTIMATED EFFECTIVE WEB THICKNESS AT THE ♦ 
C '* TERMINATION PORT IN INCHES * 


C ******** ************************************************************* 

IF (Y.LE.O.O) WRITE ( 6 » 606 ) L T F , DT P , THET TP , TAUEFF 
THE TTP= THETTP/ 5 7. 29 578 

DA BT=NT* 3. 14159*1 ( DTP+2 . *Y ) * l L TP-Y/ S I N IT HETT P ) ) - l DTP+2 . *Y )**2 /4. ♦ 
lIY*DTP/2. )*< OTP/2. )*( l.-l./S INI THETTP) )) 

IFIY.GE. TAUEFF) DABT=0.0 
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371 IF(Y.GT.O.O) GO TO 52 
IF (NT. NE. 0.0) GO TO 45 
LTP=0.0 
0TP=0.0 

45 IF(GRAIN.N£.2) GO TO 49 
IGCI=0.0 
LGNI=0.0 
01 SQO-O. 0 
00SQD=4.*RCSQD 
49 IF (GRAIN. EU. 1 ) LGSI=0.0 

VC 1= 1 . 1* ( ANUM*0 1 SUD*( LGC 1 *LGNI ) ♦ ( ANUM*OOSQO-AGSI ♦ 
1 LGS l+NT*LTP* ANUM^DTP*DTP )♦ VC IT 
52 B3P=0.0 
BBS-0.0 
BBN-0.0 

ABPORT=ABPT +AB PC*ABP S+DABT+BBP 
IF (K.LT.2) XETH=A8PC/(ABPORT*1.E-20) 
A6SL0T=ABST+AB SC+ABSS+BBS 
A8N0Z=ABNT ♦ABNC+ABNS+68N 
AB TT= ABPT+ABST+A6NT 
IF (K.GE.2 ) GO TO 55555 
SUMAB=ABPORT+ABSLQT*ABNOZ 
55555 CONTINUE 

IF(K.EQ.O) GO TO 99 

IF (K.EO. I ) ABMAIN=ABPORT*ABSLOT +ABNOZ— ABTT 
K=K*l 

IF (K.GT .2 ) GO TO 69 
GO TO 2 

69 ABTO= A8PURT ♦ABSLQ T ♦ ABNOZ-ABT T 
99 CONTINUE 

IF(Y.GT.O.O) GO TO 70 
A8P1= ABPORT 
ABN1=ABN0Z 
ABSl=ABStOT 

70 ABP2=(ABPl«-ABP0RT)/2. 

ABN2=( ABNUABNOZ )/2. 

A3S2=(ABSI+ABSL0T)/2. 

IF l INPUT. EU. I ) GO TO 76 
GO TO ( Tl» 72, 73, 74), ORDER 

71 APHEAD=APHS1 
APNOZ=APNT 
SG=SGH 

GO TO 75 

72 AP HEAO = APHT 1 
APNUZ = AP NT 
SG =0. 0 

IFIGRAIN.EQ.3) SG= ( SGH+SGN )/ 2. 

GO TO 75 
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73 APHEA0=APHT1 
APNOZ-APNS 
SG=SGN 

GO TO 75 

74 APHEA0=APHS1 
APNOZ=APNS 
SG=SGN 

GO TO 75 

76 APHEAD= APHT 
APNOZ=APNT 

75 Y= YB 

DIFF=SUMAB-SUM2 

DADY=DIFF/DELY 

ABP1=ABPQRT 

ABNl=ABNOZ 

ABSl=ABSLOT 

IF TZW.GE.O.O) GO TO 77 

ABM1=ABMAIN 

ABMAI N=ABTO 

AB TO=ABMl 

77 RETURN 

500 FORMAT!9X,I2,9X,I2,8X,I2,6X,F4.0,9X, 12, 7X, 12) 

607 FORMATI//.20X, 19HGRAIN CONFIGURATION) 

600 FORMAT! 13X,7HINPUT= , 1 2 , / , 13X, 7HGRAIN* , 12,/ , 13X, 6HSTAR= ,I2,/,I3X 
1,4HNT= «F4.J,/,13X, 7HORDER= , 12 ,/ , 13X , 5HCOP= ,12,//) 

507 FORMAT! 6X, F6. 2, 10X,E 1 1.4, 10X , E l 1 . 4, 8X, E l 1 .4, / , 22X, El l .4, 

I9X,EI1.4,8X,EII.4) 

610 FORMAT l / 13X, 40HTABULAR VALUES FOR YT EQUAL ZERO READ IN) 

583 FORMAT ( L3X, 5HA8PK=» I PE 11. 4, 5X, 5HABSK=, IP El 1.4, 5X, 5HABNK=, IPE i 1.4, 

1 5X,5HAPHK=, 1PE11.4.5X, 5HAPNK= , 1 PE 1 1 . 4 , // ) 

584 FORMAT l 1 3X , 5HVC l T= , IPE 1 1.4,//) 

505 FORMAT !6X,F 7.3 ,9X, El l .4, luX, El l .4, 8X , El l .4,/ , 22X, El 1 .4, 9X, El 1 ,4) 

611 FORMAT I/13X, 23HTABULAR VALUES FOR YT= ,F7.3,SH READ IN) 

501 FORMAT (5X,F 0.2, 6X,F7.3,9X, F7.3, 5X,F6. 2 ,9X,F8 . 5, /, 7X,F8.2, 7X,F7.2,9 
IX, F8.5,9X,F8.5) 

601 FORMAT (20X, 19HC.P. GRAIN GEOME TRY , / , 13 X , 4HD0= ,F8.2»/» 13X*4HDI= ,F 
17.3,/»13X,7HDEL0I= ,F 7. 3 , / , 1 3X , 3HS= ,Ft».2,/,l3X, 8HT HET AG= ,F9.5,/» 
2 1 3 X, 6HLGC 1= »F8.2»/»13X»6HLGNI= , F 7 . 2 , / , 1 3X , 8HTHE TON- ,F9.5,/,13X, 
38H THE TCH= ,F9.5 ,//) 

502 FORMAT! 5X.F6.2, 7X, F 8. 2 , 5X, F4. 0 , 5X ,F d. 3 , 9 X ,F 7 . 3 , 5X ,F4. 0 ) 

602 FORMAT ! 15X, 19HBASIC STAR GEOMETRY,/, 13X,4HNS= ,F6.2»/» 13X, 6HLGS 1= 
1,F8.2,/,13X,4HNP= ,F5. 0 , / , 13X , 4HRC= ,F0. 3, /, 13X, 6HF ILL = ,F7.3,/,13 
2X,4HNN= , F4.0,//) 

421 FORMAT 13!6X,F5.2) ,21 10X,F7 .5 ) ,6X,F5.2 ) 

4 ?c FORMAT 120X,20HWAG0N WHEEL GEOMETRY,/, 1 3X , 7HT AUWW= ,F6.2,/,13X, 

1 4HL 1= ,F6.2,/,13X,4HL2= ,F6. 2 , / , 1 3X, 8HALPHA1= ,F9.5,/,13X, 

2 8HALPHA2= ,F 9. 5 , / , 13X, 4HHW= ,F6.2,//) 

503 FORMAT (5X,F7.3, 7X,F7.3) 
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603 FORMAT! 20X f 23HTRUNCATED STAR GEOMETRY,/, 13X,AHRP* ,F7.3 ,/ , 13X ,6HT A 
IUS* ,F7.3,//) 

504 FORMAT !9X,F8.5,9X,F8*4«8X,F7»3) 

604 FORMAT (20X, 22 H STANDARD STAR GEOMETRY ,/, 1 3X, 8HTHETAF= ,F9. 5,/, 13X,8 
1HTHETAP= ,F9.4,/,l3X,7HTAUwS= ,F?.3,//J 

506 FORMAT(7X,F7.2,7X,F6.2,10X,F8.5,iOX,F7.3) 

606 FORMAT { 20X, 2 5H TERM (NAT ION PORT GEOMETRY, /, 13X,5HLTP= , F7.2, /, 13X, 5 
IH0TP= ,F6.2,/,13X, EHTHETTP® ,F8. 5, / , 1 3X, 6HTAUEFF= ,F 7.3,//) 

END 
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C 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE OUTPUT 

********************************************************************* 


* SUBROUTINE OUTPUT CALCULATES BASIC PERFORMANCE PARAMETERS * 

* AND PRINTS THEM OUT AS A FUNCTICN OF DISTANCE BURNED * 

* (WEIGHT CALCULATIONS ARE PERFORMED IN THE MAIN PROGRAM) * 

* T IS THE TIME IN SECS * 

* Y IS THE DISTANCE BURNED IN INCHES * 

* RNOZ IS THE NOZZLE END BURNING RATE IN INCHES/SEC * 

* RHEAD IS THE HEAD ENO BURNING RATE IN INCHES/StC * 

* PONOZ IS THE STAGNATION PRESSURE AT THE NOZZLE END IN PSIA ♦ 

* PHEAD IS THE PRESSURE AT THE HEAO END OF THE GRAIN IN PSIA * 

* PTAR IS THE PORT TO THROAT AREA RATIO * 

* MNOZ IS THE MACH NUMBER AT THE NOZZLE END OF THE GRAIN * 

* SUMAB IS THE TOTAL BURNING AREA OF PROPELLANT IN IN**2 * 

* SG IS THE BURNING PERIMETER IN INCHES OF THE STAR SEGMENT * 

* (IF ANY) * 

* PATM IS THE ATMOSPHERIC PRESSURE AT ALTITUDE IN PSIA * 

* CF VAC IS THE THEORETICAL VACUUM THRUST COEFFICIENT * 

* F VAC IS THE VACUUM THRUST IN LBS * 

* F IS THE THRUST IN LBS AT AMBIENT PRESSURE * 

* ISP IS THE DELIVERED SPECIFIC IMPULSE IN SEC AT AMBIENT * 

* PRESSURE * 

* CF IS THE THEORETICAL THRUST COEFFICIENT AT AMBIENT PRESSURE * 

* VC IS THE VOLUME OF CHAMBER GASES IN IN**3 * 

* MOOT IS THE WEIGHT FLOWRATE IN LB/SEC * 

* CFVD IS THE DELIVERED VACUUM THRUST COEFFICIENT * 

* ITOT IS THE ACCUMULATED IMPULSE IN LB-SEC OVER THE * 

* TRAJECTORY * 

* ITVAC IS THE ACCUMULATED VACUUM IMPULSE IN LB-SEC * 

* I SP VAC IS THE DELIVERED VACUUM SPECIFIC iMPULSt IN SEC ♦ 

1000 CONTINUE 

* WP IS THE EXPENDED PROPELLANT WEIGHT IN LB * 

* RADER IS THE NOZZLE THROAT EROSION RATE IN IN/SEC ♦ 

* EPS IS THE NOZZLE EXPANSION RATIO * 

* ALT IS THE ALTITUDE IN FT * 

* DT IS THE NUZZLE THROAT DIAMETER IN IN * 

* APHEAD IS THE HEAD END PORT AREA. IN IN**2 * 

* APNOZ IS THE NOZZLE ENO PORT AREA IN IN**2 * 

* COF IS THE CHARACTERISTIC THRUST COEFFICIENT * 

* CFO IS THE DELIVERED THRUST COEFFICIENT AT AMBIENT PRESSURE * 


K C * ETHETA IS THE TANGENTIAL STRAIN OF THE C.P. GRAIN AT THE BORE * 


* RH05 IS THE DENSITY OF THE PROPELLANT IN SLUGS/IN**3 AT THE * 

* CHAMBER PRESSURE AND TEMPERATURE TGR * 

* RAT P IS THE RATIO OF EXTERNAL TO INTERNAL PRESSURES UN THE * 

►C * C.P. GRAIN * 

►C * RATR2 IS THE SQUARE OF THE RATIO OF BCRE RADIUS TO OUTSIDE * 

► C * RADI US OF THE C.P. GRAIN * 

►C * XETH IS THE FRACTION OF THE TOTAL BURNING BURE SURFACE * 
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► C * ASSOCIATED WITH THE C.P. GRAIN * 

► C * YETA IS THE DISTANCE BURNED IN INCHES FOR THE ENTIRE STAR * 

►C * GRAIN AND THE C.P. GRAIN ENDS WHEN THE GRAIN DEFORMATION * 

► C ♦ MODIFICATION IS IN USE * 

C a******************************************************************** 

REAL MGEN,MDlS,MNOZ,MNl, JROCK,N,L,MEl, ME, ISP * I TOT t MU* MASS* ISP VAC 
REAL M2# MOBAR, ISP2, IT VAC, MOOT, ISPV 
COMMON/CONS T 1/ ZW, AE« AT, THE TA.ALF AN 

COMMON/CGNST2/CAPGAM,ME,BOTE ,ZETAF, TB,HB, GAME, CGAME, TOPE, ZAPE 
COMMON/ VAR I Al/Y,T, DEL Y, DELTA T , PONOZ , PHEAD ,RNOZ , RHEAD , SUMAB , PHMAX 
COMMON/ VARI A2/ ABPQRT , ABSLOT , ABNQZ, APHE AD, APNOZ ♦ DACY , ABP2 , ABN2 , ABS2 
COMMON/VARIA3/I TOT, ITVAC, JROCK, ISP, ISPVAC,MDIS,MNOZ,SG,SUMMT 
COMMON/ VARI A5/ABMAIN»ABT0,SUM0Y ,VCI , ABTT , PTR AN 

COMMON/ VAR I A6/WP2.CF.WP, RADER, EPS, VC, FLAST,TL AST, DT,PONTOT,WPi 
COMMON/ VAR I A7/T I ME, FV, ISPV, NX 

► COMMON/ VARI A9/ETHETA, RHO 5 ,RATP , PMOO , CMOD ,PMU »CMU» ALPTS ,RATR2 

► COMMON/ V AR 1 20/ YETA , XE TH, ISO 

COMMON/ I GN1/KA.K3, Uf S , RHO, L , PM I G, T 1 1,TI2 ,CSIG,U1,N1«Q2,N2 
COMMJN/PLOTT/NUMPLT! 16) , IPO,KDUM,NP, I OP 

01 ME NS I ON TPLOT 12 00) ,PNPLOT (200) ,PHPLOT( 200 ) ,FPL0T(200) ,F VPLOT 1200 
I ) , RNPLOT (200) ,RHPLOTI200) ,YBPL0T1200) , ABPLOT (200 ) .SGPLOT (200) , VCPL 
20T (200 ) 

DATA G/32.1725/ 

IF (NDUM.EQ. I ) GO TO 2 

ME 1 = 7.0 

NP =NP*l 

Y0=Y 

VC X=VC 

IF(Y.LE.O.O) M2=MD I S 
MPMAR=(M2 + Ml)IS)/2. 

SUMMT=SUMMT+MDBAK*OELTAT 
WP 1 = G*SU iMT 
WP 2= RHO* ( VC- VC I ) *G 
WP=( WPUWP2J/2. 

► IF(ISO.Eg.l) WP=WPI 
PT AR= 1 • / JROCK 

17 ME=SQRT(2./QQTE*( TOPE/2.*! AE*ME l/AT ) ** ( l ./ZAPE )- 1. ) i 
IF IABS! ME-MEI) .Lfc. 0.002) GO TO 9 
ME 1=ME 
GO TO 17 
9 CONTINUE 

PRES=(l.*80TE/2.*ME*ME)**(— GAME/BOTE) 

ALT=HB*(T/TB)**(7./3.) 

PA TM= 14.696/EXP (0. 4 J103E-04* ALT ) 

IF (MOIS.LE.O.O.OR.PON07.LE .0.0) GO TO 45 
COF=CGAME*SQRT< 2.*GAME/B0TE*i l.-PRES** (BCTE/GAME)) ) 

CF=COF*AE/AT*( PRES-PATM/PONOZ ) 

CFVAC=CF+AE/AT*PAT M/PUNUZ 
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CFO=(COF*(L.+COSIALFAN) )/2.»EPS*PRES)*ZETAF-EPS*PATM/PQN0Z 

CFVD=CFO+EPS*PATM/PONOZ 

F*COS ( THETA ) *PONOZ*AT*CFD 

IFIF.LE.O.O) F=0.0 

IF < Y.LE.0.0) F 2=F 

F8AR=IF*F2l/2. 

F VAC® COS( THETA ) *PONO/*AT*CFVD 
IF (Y.LE.0.0) FV2=F VAC 
FVBAR*(FV2*FVAC)/2. 

M0QT=MDl S*G 

I SP=F/MDQT 

ISPVAC=FVAC/MDOT 

ITOT= ITOT+F0AR*DELTAT 

ITVAC=ITVAC+FV8AR*0ELTAT 

IF(Y.LE.O.O)PON2=PONOZ 

PONBAR=( PQN2+PONOZ 1/2. 

PQNTOT®PONTOT+ POND AR*OEL TAT 

°0N2=P0N0Z 

M2®MDIS 

F2=F 

FV2=F VAC 

IF(PHEAO.GT.PHMAX) PHMAX=PHEAD 
GO TO 47 
45 CFVAC=0.0 
FVAC=0.0 
F= Q. 0 

47 WRITE (6,1) T,YB»RNOZ *RHE AD »PCNOZ »PHE AO ,PTAR, MNCZ » SUMAB » SG »PATM, CFV 
1AC»FVAC»F» I SP * CF, VCX , MOOT , CFVO , I TOT , I TV AC , I SPVAC , WP , RADER ,EP$, ALT 
2,0 T, APHEAO, APNQZ,CQF,CFD 
IF ( ISO.EQ.O) GO TO 200 

WRITE (6, 5) E THETA, RH05 » R AT P, RATR2 *XETH , Y ETA 
200 CONTINUE 

IFUPO.EQ.O) RETURN 
TPLOT (NP ) = T 
PNPLOT (NP)=PQNQZ 
PHPLOT(NP)=PHEAO 
FP LOT ( NP )=F 
FVPLOT (NP)=FVAC 
RNPLOT<NP) = RNUZ 
RHPLOT(NP)=RHEAO 
YB PLOT { NP } =YB 
AB PLOT ( NP ) = SUMAB 
SGPLOT (NP)=SG 
VC PLOT ( NP) = VC 
RE TURN 
2 NP=NP*2 
I0P=1 

DO 1004 1=1,16 
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IFINUMPLTt D.EQ.l) GO TO 1003 
GO TO 1004 

1003 GO TO (10, 20, 30, 40, 50, 55, 60, 70, 75, 80, 90, 95, 97, 100, 110, 115), I 

10 CALL PLOTITITPIQT, ‘TIME ( SECS) * * 1 l,PHPLCT , *PHEAD (PSIA)*, 12, 

1 PNPLOT » ' PONOZ* *5»NP,1, * DUMMY • ,5) 

GO TO 1004 

20 CALL PLOT IT (TPLQT , *T IME (SECS) * « 11, PNPLOT, 'PONOZ (PSIA) ' , 12,PHPL0T 
1 # • PHEAD (PSIA) * ,12,NP,1, *OUMMY* , 5) 

GO TO 1004 

30 CALL PLOTlT(TPlOT,*TIME ( SECS ) • » 11 , PHPLCT, • PHEAD* , 5, PNPLOT 
1, ‘PONOZ* ,5, NP, 3, 'PRESSURE (PSIA)’ ,15) 

GO TO 1004 

40 CALL PLOTITITPLOT, ’TIME i SECS ) ’ , l l ,RHPLOT, • RHEAD (IN PER SEC)’, 18, 
1PHPL0T,* PHEAD (PS I A) • , 12 ,NP, 1, ’DUMMY’, 5) 

GO TO 1004 

50 CALL PLOTIT( TPLOT , ’TIME ( SECS) • , l 1 , RNPLOT , • RNOZ (IN PER SEC)*, 17, 

1 PNPLOT,' PONOZ (PSIA)’, 12, NP,1,’ DUMMY ’,5) 

GO TO 1004 

55 CALL PLOTlT(TPLOT,'TlME ( S EC S > • , 1 1 , RHPIOT , • RHEAD ' , 5, RNPLOT , 

1 • RNOZ* ,4, NP, 3, ’BURNING RATE (IN PER SEC)’, 25) 

GO TO 1004 

60 CALL PLOTITITPLOT, 'TIME (SECS )', 11 , ABPLCT, ’TOTAL BURNING AREA (SQ 
1 IN) ’,26, PNPLOT, 'PONOZ* ,5, NP, 1, • DUMMY* ,5) 

GO TO 1004 

70 CALL PLOTITITPLOT, 'TIME ( SEC S )*, 1 1 , SGPLOT ,’ STAR PERIMETER (IN)’, 19 
1, PNPLOT, ’ PONOZ* ,5,NP, l, 'DUMMY* ,5) 

GO TO 1004 

75 CALL PLOTIT(TPLOT, ’TIME ( SEC S ) • , 1 1 , ABP LC T , • TOT AL BURNING AREA (SQ 
1IN)*, 26, SGPLOT, 'STAR PERIMETER ( IN ) • , 19, NP, 2 , ' DUMMY • , 5 ) 

GO TO 1004 

80 CALL PLOT IT (TPLOT , 'TIME 1 S ECS )*, 1 1 , F PLOT ,’ THRUST (LBS) ’, 12, PNPLOT , 
1 • PONOZ ', 5, NP,1,’ DUMMY* ,5) 

GO TO 1004 

90 CALL PLOTITITPLOT, ’TIME ( SECS )’, 1 1 , FVPLO T ,’ VACUUM THRUST (LBS)*, 19 
1, PNPLOT, 'PONOZ* ,5,NP, 1,’ DUMMY’ ,5) 

GO TO 1004 

95 CALL PLOTITITPLOT, 'TIME ( SECS )’, 1 1 , F PLOT ,' THRUST ’, 6 , FVPLOT , 

1 ’VACUUM THRUST’ ,13, NP, 3, 'THRUST (LBS>’,12) 

GO TO 1004 

97 CALL PLOTIT ITPLUT, ’TIME ( SECS) ’, l 1 , VCPLO T ,’ CHAMBER VOLUME (IN**3)* 
1 , 22 , PNPLOT , • PONOZ • , 5 , NP , 1 , « DUMMY • , 5 ) 

GO TO 1004 

100 CALL PLOTIT IYBPLOT , 'BURNED DISTANCE ( I N ) • , 20, ABPLOT, • TOTAL BURNING 
1 AREA I SQ IN)’, 26, PNPLOT, ’PONOZ’, 5, NP,1, ’DUMMY’, 5) 

GO TO 1004 

110 CALL PLOUTIYBPLOT , 'BURNED DISTANCE ( IN) *, 2(5, SGPLOT, ’STAR PERIMETE 
IR I IN)’, 19, PNPLOT, 'PONOZ', 5, NP , 1 , • DUMMY’ ,5) 

GO TO 1004 


- 111 - 



Table C-3 (Cont'd) 


115 CALL PLOTIHYBPLOT , 'BURNED DISTANCE UN) • ,20. A8PLOT, • TOTAL BURNING 
1 AREA (SQ IN)' ,26,SGPLCIT, 'STAR PERIMETER ( IN ) • , 19.NP ,2 , ' DUMMY* , 5) 
1004 CONTINUE 
RETURN 

l F0RMATI13X.6HTIME* ,F7.2 , 12X ,3HY* ,F6.2, /, 13X.6HRN0Z* , IPE 1 1 » 4. 9H 
l RHEAO* .IPEll. 4, 9H PONOZ* .IPEll. 4, 9H PrlEAD* , 1PE ll .4, /. 13X ,6HP 
2TAR* , 1PE1 1 .4, 9H MNOZ* .IPEll. 4, 9H SUMAB* , IPEll. 4.9H $G* , 

31PE11.4./.13X.6HPATM* .IPEll. 4, 9H CFVAC* .IPEIU4.9H FVAC* ,IPE 
411.4, 9H F* »lPEll.4j/»l3X,6H ISP* .IPE11.4.9H CF* .IPEll. 

54, 9H VC* , IPEll. 4,9H MCOT* , 1PE 1 1. 4,/, 13X.6HCFVD* .IPEll. 4, 9 
6H ITOT* , IPE 11.4, 9H ITVAC* .IPEll. 4, 9H ISPVAC* , l PE 1 1 . 4, / , 13X , 6 
7HWP* .IPEll. 4, 9H RADER* .IPEll. 4, 9H EPS* .IPEll. 4, 9H ALT* 
8 , !PE11.4,/,13X,6HDT* ,iPEli.4,9H APHE AD= .IPEll. 4, 9H APNOZ* ,1 
9PE11.4.9H COF* , IPE l l .4, / , 1 3X , 6H CFO* .IPEll. 4,/) 

5 FORMAT! 13X,8HETHETA= .IPEll. 4, 9H RH05* , IPEll. 4, 9H RATP* ,IPE1 
11.4.9H RATR2* .IPEll. 4, 9H XETH* .IPEll. 4, 7H YETA* .IPEll. 4,//) 
END 
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SUBROUTINE IGNITN 

C ***** ********** <*** ******** ************** ***************************** 

C * SUBROUTINE IGNITN CALCULATES THE PRESSURE RISE OURING * 

C * THE IGNITION PERIOD * 

C * ASIG IS THE IGNITER THROAT AREA IN IN**2 * 

C ♦ WIGTOT IS THE TOTAL WEIGHT OF THE IGNITER PF.OPELL ANT IN LBS * 

C * MIGAV IS THE IGNITER AVERAGE MASS FLOW RATE OVER THE FIRST * 

C * HALF OF THE IGNITER BURNING TIRE IN LaS/SEC * 

C ♦ PC I G IS THE IGNITER PRESSURE IN LBS/IN**2 * 

C ********************************************************************* 
REAL K(4) ,L,KA,KB, JROCK, J2 , M I G, MI GA V , MSR M , ME , MDI S , MNOZ , MNO Z I , MN1 
REAL N1 , N2 , M (GAVE 

COMMON/CONS T 1/ZW,AE, AT, T HE TA,ALF AN 

COMMON/C ONST 2/ CAP GAM , ME, BOTE ,ZETAF,TB,HB, GAME, C GAME »T0PE,ZAPE 
COMMON/ VAR I A1/Y,T I G, DEL Y , DEL TAT , PCNOZ , PH EAD , RNOZ , RHE AD , SUMAB , PHMAX 
COMMON/ VAR I A2/ A6PORT , ABSLOT , ABNQZ , APHE AD , APNOZ, DAOY , A9P2, ABN2,AES2 
COMMUN/VARI A3/ I TOT , I T VAC , J ROCK , ISP, I SP V AC , MO I S , MNOZ , SG , SUMMT 
COMMON/ VARI A5/ AoMA I N , AB TO , SUMOY , VC I , AB TT , PT R AN 
COMMON/ IGN 1 /KA,Kb, UFS, R HO,L, PMIG,TI 1,T 12 , CS I G , Q L , N 1 , Q2 , N2 
COMMON/ I GN2/ ALPHA,BETA,PBIG,RRIG,OELTIG,X,TOP,ZAP 
COMMON/PLOTT/NUMPL T ( 16) , IPO.NDUM, IPT, I CP 
DIMENSION B ( 9 » 

DATA Al, A2, A3, A4/. 17476 ,-.55 1481, 1.205 536,. 171 185/ 

OAT A BID, B(2), 8(31,8(41,6151/0. ,.4, .455737,1. ,.296978/ 

DATA B(6 ) »6( 7) ,3(8 ) , ft I 9 ) / . 1 58 76 , . 2 1 o 1 , -3 .050965 , 3 . 832 864/ 

C ****♦****** + *♦****♦♦****#******#*********♦*♦♦*#**<1******************* 
C * THE A 1 S ANU B‘S ARE CONSTANTS FOR THE RUNGE-KUTTA INTEGRATION * 
C ********************************************************************* 
OATA G/3 2.172 5/ 

XXX=.05*P0N0Z 
IP LUG= 0 
PONOZI- ' '‘NGZ 
RHEADJ ‘ : AD 

RN CZ I = >L 

PHEADI^RHEAD 
DE LT T = OEl TAT 
DI SM-MD I S 
DELTAT^DELT tG 
SUMAB l=SUMAB 
MNCZl =MNOZ 
MNQZ=0.0 
RH FAD=0. 0 
RN0Z=0.0 
MO I $= 0.0 
AB 1=0,0 
T I GI =0. 0 
PC 1=14.696 
T I G = 0 ,0 
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PC NEW* 14.696 

SUMAB=0.0 

PC IG= 14.696 

PHEA0*14.696 

PONOZ* 14.696 

$L0PE*SUMA8l/L 

G2*CAPGAM*CAPGAM 

J2*JRQCK*JR0CK 

GJ*G2*J2/2. 

MIGAV*.2*AT/G 

AS IG*4.*MIGAV*CSIG/I4«*PMIG“RRIG*ITI2— TI 1) ) 

MI GTOT*G*MI GAV4|!»*PC TI2— TI l)/6. ) 

MIGAVE*NIGAV*G 

WRITE (6*999) AS IG« Ml GTOT ,M IGAVE 
WRITE! 6« 10) 

18 NNN*0 

WRITE (6« 30) PC IG 
CALL OUTPUT 
9 CONTINUE 

00 8 N*l,4 

IF IN.EQ. 1 I PC«PCI 
IFIN.EQ.2) PC«PCI+B(2)*Ktl) 

IF IN.EQ. 3) PC* PC I + 8 ( 5)*KI1)*8I6)*KI2) 

IF IN.EQ. 4) PC*PCH-8I7)*Kt 1 H-Bt 8 )*K(2)*B( 9) *K I 3) 

TIG=TIG1 *Bi N)*DELT IG 
SUMA8=AB l+SLOPE*UF S*81N) *OELT IG 
IF (SUHA8.GT . SUMA8I ) SUMAB=SUMA9l 
PHEAO=PC 

IFIMOIS.NE.O.O) PH6AD=PC*I l.+GJ) 

IF IPHEAO.LE .PTRAN) RHEA0=Q1*PHEA0**N1 
IF IPHEAD.GT.PT RAN )RHEA0=Q2*PHEAD**N2 
IF I T IG.LE.T 1 1 ) PC I G=PM IG*T IG/T 1 1 

IFCTIG.GT.T I1.AN0.PCIG.GT.PH6A0J PC I G=PM IG-KR IG* I T IG-T 1 1 ) 
IFIPCIG.LE.PHEAD) PCIG=PHEAD 
MI G=0.0 

IF IPCIG.GT.PHEAD.AND.T IG.LE.T 12/2.) MIG=PCIG*ASIG/CS IG 

CSTR=KA+K8*PC 

MD I S= PC* AT /C ST R 

IFIPC.LE.PBIG.AND. IPLUG.EO.O) GO TO 7 

1 PLUG* l 
PNUZ=MNOZI 
PNOZ=PC*II.-GJ) 

ZI T=MDI S*X/APNOZ 
RNl=KHEAO 
AZ=ALPHA*ZIT**.8 
XL=UFS*TIG 
IFIXL.GT.L) XL = L 

4 EX = XL**.2*EXPl BETA*RN1*PHU/Z IT) 
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IFIPNOZ.LE.PTRAN|RNOZ^RN1-(RN1-Q1*PNOZ**N1-AZ/EXI/U.4-AZ*BETA*RHO/ 
2(1 1 T*EX) ) 

IF ( PNOZ.GT .PTRAN)RN0Z=RN1— IRNi-Q2*PNOZ**N2-AZ/EX) /( l.*AZ*BE TA*RHO/ 
2IZ IT*EX) I 

IF (ABS( RN1—RN0Z)-LE. 0.002) GO TO 5 

RNl=RNOZ 

GO TO 4 

7 MDIS=0.0 
HNQZ=0.0 
PNOl-PC 
RNOZ=RHEAO 

5 CONTINUE 

MSRM=RHO*SUMAB*lRNOZ*RHEADI /2. 

OE NOM=tVCI/U2.*CSTR*C$TR*G2))*Ci.-{2.*Ke*PC)/CSTR) 

OPDT= ( MIG+MSRM— MO I SI /OENQM 
IFIOPOT.LT. 0.0. aNO.PC.LT. 20.0) DP0T-=0 .0 
KI N)=0£LT I G*DPDT 

8 CONTINUE 

PCNEW=PCI*A1*K< 1 )*A2*K(2)*A3*KI3)*A4*K(4) 

PHEAO=PCNEW 

IF (MDI S.GT .0.0 ) PHEAD=PCNEW* ( L.+GJ) 

PQNOZ=PCNEW 

XX Y-A3SI PONOZ-PONOZ I ) 

IF( PC NEW. LE. 1.001* PC I. AND. SUMAB. EU.SUMAB I • ANQ.XXY *LE .XXX) GO TO 13 

AB 1 = SUMAB 

TIGI=TIG 

PC I=PCNE W 

NNN=NNN«-1 

IF (NNN.GE.5 I GO TO 18 
GO TO 9 
13 CONTINUE 
CALL OUTPUT 
WR ITE ( 6 » 30 ) PC I G 
OELTAT=DELTT 
MO IS=0 I SM 
SUMAB = SUMAB I 
PONOZ=PGNUZI 
RHfcAO=RHEAO I 
P.NQZ=RNUZ I 
PHEAD=PHEAOI 
MNUZ=MNUZl 

IF (IPO.NE.2.ANO.IPO.NE.3I GO TO 53 
NOUM= 1 
CALL OUTPUT 
N0UM=0 
53 CONTINUE 
I P T=0 
RETURN 
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999 FORMAT!/// «20X » 25 H IGNITER SIZE CALCUL AT I CNS. / . 13X , 5HAS IG= . F 7 . 2t /. 
1 13X»7HMIGT0T=*F7.2»/»13X»6HMI GAV= »F8.3»///) 

10 FORMATI33X»2 8H ******************* **»*♦♦♦ ***/»33X» 28H**** IGNITION 
l TRANS I ENT **** , /*33X,28H* *************************** | 

30 FORMAT ( 1 3X» 6HPCIG= .IPEU.4I 
ENO 


SUBROUTINE INTRP 1 ( Y , T ,N , TT ,DY , ICHK » 

DIMENSION YIN) , TIN) 

N1=N— 1 
07=0.0 

IF ( ICHK) 2*2*3 

2 DO 1 1=1 »N1 

IFITT.GE.TI U.AND.TT.LT.TlIill) DY=(( Y l I ♦ 1 )- YC III /I TCI ♦ i ) — T ( 1 ) ) ) 
2*1 TT-TU ) )*YI I ) 

IFIOY.NE.O.O) RETURN 
1 CONTINUE 

3 DO 4 l =1 *N1 

IF4TT.LE.T (I ) . AND. TT.GT.T I 1 + 1) ) 0Y=( (YI i ♦ l )- Y( I ) ) /(T( I +1 )~Tl ()) ) 
2*( TT-T4 I )UY{ I) 

IF { OY.NE.O.O) RETURN 

4 CONTINUE 
RETURN 
END 
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SUBROUTINE PLOT IT ( X.XHOR.KX, Y, YHDR, NY, T. THOR , NT ,NP .NPLOT , DUMMY , NO ) 
******** ************************************************************ 

* SUBROUTINE PLOT IT PLOTS TWC DEPENDENT VARIABLES* Y AND T. 

* VERSUS AN INDEPENDENT VARIABLE. X 

* XHDR* YHDR* ANO THDR ARE THE HEADINGS FOR THf X, Y, AND T 

* AXES. RESPECTIVELY 

* KX, NY. ANO NT ARE THE NUMBER OF CHARACTERS IN THE X, Y, AND 

* T AXES HEADINGS, RESPECTIVELY (MAX OF 32 IN EACH) 

* NP IS THE NUMBER OF POINTS TO BE PLCTTEO PLUS 2 

* VALUES FOR NPLOT ARE 

* 1 FOR Y ONLY PLOTTEO VERSUS X 

* 2 FOR Y AND T PLCTTEO VERSUS X ON SAME AXES 

* WITH INDIVIDUAL SCALES 

* 3 FOR Y AND T PLCTTEO VERSUS X ON SAME AXES 

* WITH SAME SCALES 

* OUMMY IS THE HEAOING FOR THE DOUBLE AXIS (NPLOT-3) 

* NO IS THE NUMBER OF CHARACTERS IN OUMMY 
******************************************************************** 

DIMENSION XHOR (8)»YHDR(8) , THDR 1 8 ) .DUMMY ( 8),X{NP),Y(NPJ,T(NP) 

NX=-K X 

NM=NP-1 

NN=NP-2 

IF(NPLOT.EQ.i) GO TO 9 
CALL SCALE! T, 4. »NN» I ) 

9 CALL SCALE ( X, 8. .NN . 1 ) 

CALL SCALE( Y,4. »NN» 1 ) 

IF (NPLOT. NE. 3) CALL AX I S ( 0. , 0. , YHDR ,NY ,4. , 180. , Y! NM| , Y (NP ) ) 

IF (NPLOT . EQ. 3 ) CALL AX 1 S ( 0 ., 0. , DUMMY .NO , 4. , 1 80. , Y ( NM ) , Y ( NP)) 

CALL AXiS(Q.,0.,XHDR,NX,8.,90.,X(NM) ,X(NP)) 

I F ( NPLOT • EQ. 1 ) GO TO 12 
DO 11 1=1, NN 

11 T ( I )=-T( I ) 

12 00 13 1=1, NN 

13 Y ( I ) =-Y ( I ) 

CALL LINE(Y.X,NN, 1,0,1) 

CALL PLOT (0..0..3) 

IF( NPLOT. EQ.l) GO TO 24 

IF (NPLOT .E0.2 ) CALL PLOT (0 . 5, 2 ) 

IF (NPLOT. EO. 2 ) CALL AX l S ( 0. ,-. 5 , THDR ,N T , 4. , 1 80. , T ( NM ) , T ( NP ) ) 

CALL L INE(T,X,NN, 1 ,0,2) 

DO 25 1=1 »NN 

25 T ( I ) = -T ( I ) 

24 DO 26 1=1, NN 

26 Y( I ) = — Y { I ) 

IF (NPLOT .EQ. I ) GO TO 32 

CALL SYMB0L(-4.35,.52, ,1 ,l,0.,0) 

CALL SYMBOL (-4.2, . 52, . 1,2, 0. ,0) 

CALL SYMBOL (-4. 3,. 65,. 1. YHDR, 90., NY) 

CALL SYMBOL (-4 . 1 5 , . 65 , . 1 , THDR , 90 . ,N T ) 

32 CALL PLOT (8. 5,0., -3) 

RETURN 

END 
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